Orbits in Gridded Potentials
Here is another example how to compute orbits
within NEMO, but now within a potential that is defined on a grid. This
opens up an way to take real data (yes yes, from the sky).
First we need a small excursion into the obervers domain. Here are the basic
steps (and you can skip steps if they've been taken already):
- import your fits data into miriad, using fits in= out= op=xyin
- make sure the center of the galaxy is at the reference pixel,otherwise
fake it by using puthd in=map/crpix1 value= and
puthd in=map/crpix2 value=. This will mess up the RA/DEC,
but we only care about relative numbers here.
- deproject the image to face on. Since the reference pixel has been
set properly in the previous step, you only need to decide on
the position angle and the inclination.
deproject in= out= pa=xxx inc=yyy
- use any smoothing or transformation you like to make a more realistic (?)
density distribution. Perhaps even using ellint.
- calculate the potential from the density and with an assumed scale height
of the gas: potfft in= out= h=zzz. The units of scale height
are the same units as the pixel, unless the distance in Mpc is given,
these need to be given in pc.
- hi saf
- output the file in fits for importing into nemo: fits in= out= op=xyout
There are three ways to make an
that can be used with the potname=ccd
You may want to test some values of this potential
or simply display the potential in ds9 with the
- Import the fits file into NEMO using
for example from a file you exported in MIRIAD
fitsccd in=pot1.fits out=pot1
- Use ccdmath
to define an image from scratch. Here is an
example of a Plummer potential:
ccdmath out=pot2 fie="-1/sqrt((%r/100)**2+1)" size=800,800,1 cdelt=0.01,0.01 crpix=400,400
Note the peculiar normalization of the radius because the pixel size
to define an image from another potential
descriptor, e.g. the same plummer potential as in the previous
potccd pot3 plummer x=-4:4:0.01 y=-4:4:0.01
Now create an orbit wit
mkorbit orb1 x=1 y=0 z=0 vx=0 vy=1 vz=0 potname=plummer
double IOM -0.207107 1.00000 0.00000
where for convenience the Plummer potential was also given (not
required if you pass all 6 phase space coordinates) to check
if the energy is bound and negative (-0.207107).
Now the orbit can be integrated using
using either potential descriptor:
orbint orb1 orb1.out 10000 0.01 potname=plummer
orbint orb1 orb2.out 10000 0.01 potname=ccd potfile=pot2
orbint orb1 orb3.out 10000 0.01 potname=ccd potfile=pot3
and plotted using
orbplot orb1.out yapp=1/xs xrange=-4:4 yrange=-4:4
orbplot orb2.out yapp=2/xs xrange=-4:4 yrange=-4:4
orbplot orb3.out yapp=3/xs xrange=-4:4 yrange=-4:4
and you should see they are basically all the same.
potcode program can compute
orbits in static potential with certain forms of interaction between