Back to the graduate course page 

Problem set #1. Due  ................. except for the Poisson solver which is due .................

(BT=Binney and Tremaine 1987). [Ignore any $ and \ signs in the text of the problems, they may be some remains of the Tex/Latex way of marking the mathematical expressions I once tried to propose, but since then I tried to change everything into ASCII again.]
Problem #1. This is a warm-up problem. The bulge of NGC1234 has a bulge potential of the form
Phi(r) = -GM/(r+r_0), where r_0 and M are constants.
Investigate thoroughly the properties of the bulge: density rho(r), mass M(r), total mass M_{tot}, W (pot. energy), v_c(r) (=circular speed, or rotation curve). Plot the functions and determine their asymptotic behavior (ie to what functions, not just values, they tend). What's the relation between M and M_{tot}. Compare rusults with the Plummer-Kuzmin model.

The rotation curve of the whole galaxy is known to rise from 220 km/s at the center to a maximum of 300 km/s at r=3 kpc, and to further decline smoothly to 220 km/s again. Determine all the relevant parameters of the bulde and halo of NGC1234 and compare the densities of its components as a function of radius. Plot the mass-to-light ratio Upsilon(r), assuming that the halo is absolutely dark and Upsilon=1 for the bulge.

Problem #2. Solve probl. BT [2-1], and [2-2] (unless you already did during the undergraduate course). 
Problem #3. Guided by chapt. 2.6.3 of BT, redo the authors' derivations, done in the cylindrical coordinates, in the much simpler Cartesian system (x,y,z). Use the following basis function (harmonic waves, or plane waves):

(some function of z)  *  exp i(k_x * x + k_y * y)
or, if you prefer,  
exp i(k_x * x + k_y * y + some funct. of z )

Your results should be of a form analogous to these on p.76 of BT, with Hanckel transforms replaced by Fourier transforms of rho(x,y)

Problem #4. Generalize the result of problem 3 from the case of a razor-thin but otherwise arbitrary density distribution Sigma(x,y), to a general case of arbitrary rho(x,y,z).

Hint: a 3-D object consists of slices of 2-D distributions of the form Sigma(x,y)/dz.

Then try to look back at ALL the methods and equations in BT for finding Phi(x,y,z) and generalize them to a finite-thickness disk case. 

Problem #5. Write a C or Fortran subroutine which given an input 2D or 3D array of density values (in rectangular coordinates) returns the gravitational potential, or separate arrays with force components Fx, Fy etc. (Why is that a preferred output if your final aim is finding force not potential?). The subroutine should use one of the fast transform techniques. It should also include a simple direct summation subroutine to check the correctness of the fast method. What speed-up factor did you achieve by using FFT instead of direct summation?

Hint: Here is the Fortran source code for the many-dimensional FFT routine from Numerical Recipes. You can highlight it and read it into your editor or dump this whole page by selecting Save As from File menu of your Netscape. Please read the appropriate chapter of the above book to find out how it works, what is the structure of input and output, and how to avoid `aliasing' while using the Fourier transforms. (You may download the relevant section of the book in Postscript or Adobe Acrobat format from the Numerical Recipes books On-Line. This routine has long lines, consult with me in case of compilation problems; also notice it is fed single-precision data!]

One more note: there is at least one subroutine that may be faster in case of real data, in the recent editions of Num. Rec. Try using it for real inputs.  

        SUBROUTINE FOURN(DATA,NN,NDIM,ISIGN)
        REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
        DIMENSION NN(NDIM),DATA(*)
        NTOT=1
        DO IDIM=1,NDIM
                NTOT=NTOT*NN(IDIM)
        ENDDO
        NPREV=1
        DO IDIM=1,NDIM
                N=NN(IDIM)
                NREM=NTOT/(N*NPREV)
                IP1=2*NPREV
                IP2=IP1*N
                IP3=IP2*NREM
                I2REV=1
                DO I2=1,IP2,IP1
                        IF(I2.LT.I2REV)THEN
                                DO I1=I2,I2+IP1-2,2
                                        DO I3=I1,IP3,IP2
                                                I3REV=I2REV+I3-I2
                                                TEMPR=DATA(I3)
                                                TEMPI=DATA(I3+1)
                                                DATA(I3)=DATA(I3REV)
                                                DATA(I3+1)=DATA(I3REV+1)
                                                DATA(I3REV)=TEMPR
                                                DATA(I3REV+1)=TEMPI
                                        ENDDO
                                ENDDO
                        ENDIF
                        IBIT=IP2/2
1                       IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN
                                I2REV=I2REV-IBIT
                                IBIT=IBIT/2
                        GO TO 1
                        ENDIF
                        I2REV=I2REV+IBIT
                ENDDO
                IFP1=IP1
2               IF(IFP1.LT.IP2)THEN
                        IFP2=2*IFP1
                        THETA=ISIGN*6.28318530717959D0/(IFP2/IP1)
                        WPR=-2.D0*DSIN(0.5D0*THETA)**2
                        WPI=DSIN(THETA)
                        WR=1.D0
                        WI=0.D0
                        DO I3=1,IFP1,IP1
                                DO I1=I3,I3+IP1-2,2
                                 DO I2=I1,IP3,IFP2
                                   K1=I2
                                   K2=K1+IFP1
                                   TEMPR=SNGL(WR)*DATA(K2)-SNGL(WI)*DATA(K2+1)
                                   TEMPI=SNGL(WR)*DATA(K2+1)+SNGL(WI)*DATA(K2)
                                   DATA(K2)=DATA(K1)-TEMPR
                                   DATA(K2+1)=DATA(K1+1)-TEMPI
                                   DATA(K1)=DATA(K1)+TEMPR
                                   DATA(K1+1)=DATA(K1+1)+TEMPI
                                 ENDDO
                                ENDDO
                                WTEMP=WR
                                WR=WR*WPR-WI*WPI+WR
                                WI=WI*WPR+WTEMP*WPI+WI
                        ENDDO
                        IFP1=IFP2
                GO TO 2
                ENDIF
                NPREV=N*NPREV
        ENDDO
        RETURN
        END
If you managed to solve all this and have nothing to do... see the next set.