https://github.com/florentrenaud/nbody6tt
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
imf.f
SUBROUTINE IMF(BODY10,BODYN)
*
*
* Initial mass function.
* ----------------------
*
INCLUDE 'common6.h'
*
*
* Generate Scalo IMF (Eggleton, Fitchett & Tout, Ap.J. 347, 998).
ITER = 1
* Set upper limit of mass spectrum & corresponding mass fraction.
BODYM = BODY10
XM = 0.998
* Find initial value of XM for upper & lower mass by iteration.
1 Y0 = (1.0 - XM)**0.75 + 0.032*(1.0 - XM)**0.25
Y2 = (0.75 + 0.008/SQRT(1.0 - XM))/
& ((1.0 - XM) + 0.032*SQRT(1.0 - XM))
Y1 = BODYM - 0.19*XM/Y0
YPR = -0.19*(1.0 + XM*Y2)/Y0
* Try next guess of Newton-Raphson iteration.
XM = XM - Y1/YPR
* Check upper & lower limits.
IF (XM.GT.1.0) XM = 0.99999
IF (XM.LT.0.0) XM = 0.001
* Set current mass and check the convergence.
BODYS = 0.19*XM/Y0
IF (ABS(BODYS - BODYM).GT.1.0D-06*BODYM) GO TO 1
*
* Save upper value of XM and perform second iteration.
IF (ITER.EQ.1) THEN
X1 = XM
BODYM = BODYN
XM = 0.4
ITER = 2
GO TO 1
END IF
*
* Set incremental fraction and assign individual masses.
DX = (X1 - XM)/FLOAT(N - 1)
ZMASS = 0.0D0
DO 10 I = 1,N
XI = X1 - DX*FLOAT(I-1)
ZM0 = (1.0 - XI)**0.75 + 0.032*(1.0 - XI)**0.25
BODY(I) = 0.19*XI/ZM0
ZMASS = ZMASS + BODY(I)
10 CONTINUE
*
WRITE (6,20) BODY(1), BODY(N), ZMASS/FLOAT(N)
20 FORMAT (/,12X,'REALISTIC MASS FUNCTION:',' BODY(1) =',1PE9.2,
& ' BODY(N) =',E9.2,' <M> =',E9.2)
*
* Replace input value by actual mean mass in solar units.
ZMBAR = ZMASS/FLOAT(N)
*
RETURN
*
END