//package:  starlab
//date:     11-feb-95
//author:   Piet Hut
//task:     mkplummer
//does:     constructs a Plummer model, with a spatial or mass cut-off
//index:    model, init
//also:     mkcube
//history:
// 
//    version 1:  July 1989   Piet Hut               email: piet@iassns.bitnet
//                            Institute for Advanced Study, Princeton, NJ, USA
//    version 2:  Dec 1992   Piet Hut  --  adopted to the new C++-based starlab
//+files:
//          mkplummer.C
//
//-files:
//+desc:
//     Mkplummer builds a nbody system according to a Plummer model, in virial
//  units (M=G=-4E=1, with E the total energy), and finite spatial extent 
//  which can be regulated by specifying mfrac or rfrac or using their default
//  values.                
//
//     litt: S.J. Aarseth, M. Henon and R. Wielen (1974),
//           Astron. and Astrophys. 37, p. 183.
//-desc:
/*
 *.............................................................................
 */

#include "dyn.h"

#ifdef TOOLBOX

#define  MFRAC_DEFAULT  0.999                 /* radial cutoff               */
#define  RFRAC_DEFAULT  22.8042468            /* corresponding radial cutoff */
                                              /*   (not used right now)      */
#define  SEED_STRING_LENGTH  60

//+note:
//       after sprinkling in particles according to a Plummer
//       distribution, the whole system is shifted in position
//       and velocity so as to put the center of mass at rest
//       at the coordinate center.
//-note:
//+bugs:
//	None.
//-bugs:
/*
 *-----------------------------------------------------------------------------
 */
local void  mkplummer(dyn * b, int n, real mfrac, real rfrac)
    {
    int  i;

    real  partmass;		   /* equal mass of each particle	     */
    real  radius;		   /* absolute value of position vector      */
    real  velocity;		   /* absolute value of velocity vector      */
    real  theta, phi;		   /* direction angles of above vectors      */
    real  x, y;		           /* for use in rejection technique         */
    real  scalefactor;             /* for converting between different units */
    real  inv_scalefactor;         /* inverse scale factor                   */
    real  sqrt_scalefactor;        /* sqare root of scale factor             */
    real  mrfrac;                  /* m( rfrac )                             */
    real  m_min, m_max;            /* mass shell limits for quiet start      */
    dyn * bi;

    b->set_mass(1.0);
    partmass = 1.0 / ((real) n);

//+desc:
//  Calculating the coordinates is easiest in STRUCTURAL units;
//  conversion to VIRIAL units will be performed below.
//
//    Recipe for scaling to the proper system of units:
//
//  Since G = M = 1, if we switch from a coordinate system with
//  length unit  r_old  to a coordinate system with length unit  r_new ,
//  the length units simply scale by a factor  C = r_new / r_old .
//  Consequently, the coordinate values of physical quantities
//  such as positions should transform inversely to maintain the same
//  coordinate-invariant meaning. Similarly, the square of the velocities
//  should transform inversely proportional to the positions,
//  since  GM = 1  (cf. a relation such as  v*v = G*M/r ).
//<pre>
//  To summarize: If
//                       r_unit(new) = C * r_unit(old)  ,
//                then
//                       pos(new) = (1/C) * pos(old)
//                and
//                       vel(new) = sqrt(C) * vel(old)  .
//</pre>
//-desc:

    scalefactor = 16.0 / (3.0 * PI);
    inv_scalefactor = 1.0 / scalefactor;
    sqrt_scalefactor = sqrt( scalefactor );

/*
 *  Now convert  rfrac  into an equivalent mfrac, if necessary:
 */
    if (rfrac < VERY_LARGE_NUMBER) { // Note: the following powers can cause
	                             // problems on some compilers...

	rfrac *= scalefactor;        // Convert from VIRIAL to STRUCTURAL units
	mrfrac = rfrac*rfrac*rfrac / pow(1.0 + rfrac*rfrac, 1.5);
	if (mrfrac < mfrac)
	    mfrac = mrfrac;          // mfrac = min(mfrac, m(rfrac))

    }

/*
 *  now we construct the individual particles:
 */
    for (i = 0, bi = b->get_oldest_daughter(); i < n;
         i++, bi = bi->get_younger_sister())
	{
	bi->set_mass(partmass);
//+desc:
//  the position coordinates are determined by inverting the cumulative
//  mass-radius relation, with the cumulative mass drawn randomly from
//  [0, mfrac]; cf. Aarseth et al. (1974), eq. (A2).
//-desc:
        m_min = (i * mfrac)/n;
        m_max = ((i+1) * mfrac)/n;
	real rrrr = randinter(m_min,m_max);
	radius = 1.0 / sqrt( pow (rrrr, -2.0/3.0) - 1.0);
	theta = acos(randinter(-1.0, 1.0));
	phi = randinter(0.0, TWO_PI);

        bi->set_pos(vector(radius * sin( theta ) * cos( phi ),
			   radius * sin( theta ) * sin( phi ),
			   radius * cos( theta )));
//+desc:
//  the velocity coordinates are determined using von Neumann's rejection
//  technique, cf. Aarseth et al. (1974), eq. (A4,5).
//  First we take initial values for x, the ratio of velocity and escape
//  velocity (q in Aarseth et al.), and y, as a trick to enter the body of the
//  while loop.
//-desc:
	x = 0.0;
	y = 0.1;
/*
 *  Then we keep spinning the random number generator until we find a pair
 *  of values (x,y), so that y < g(x) = x*x*pow( 1.0 - x*x, 3.5) . Whenever
 *  an y-value lies above the g(x) curve, the (x,y) pair is discarded, and
 *  a new pair is selected. The value 0.1 is chosen as a good upper limit for
 *  g(x) in [0,1] : 0.1 > max g(x) = 0.092 for 0 < x < 1.
 */
	while (y > x*x*pow( 1.0 - x*x, 3.5))
	    {
	    x = randinter(0.0,1.0);
	    y = randinter(0.0,0.1);
	    }
/*
 *  If y < g(x), proceed to calculate the velocity components:
 */
	velocity = x * sqrt(2.0) * pow( 1.0 + radius*radius, -0.25);
	theta = acos(randinter(-1.0, 1.0));
	phi = randinter(0.0,TWO_PI);

        bi->set_vel(vector(velocity * sin( theta ) * cos( phi ),
			   velocity * sin( theta ) * sin( phi ),
			   velocity * cos( theta )));
	}

//  Transform to center-of-mass coordinates, and scale to standard parameters,
//  assuming a softening parameter of zero.

    b->to_com();

    if (n > 1) {

        real kinetic, potential;

	get_energies(b, 0.0, kinetic, potential);
	scale_virial(b, -0.5, kinetic, potential);
	scale_energy(b, -0.25, kinetic + potential);
        }
    }


//usage:    mkplummer -n # [options]
//parameters:
//  By default input is read from standard input, output to standard output.
//
//option:   -n#
//                   where # is the number of bodies in the Nbody system.
//                   There is no default, this is a required parameter.
//option:   -s # (seed)
//                   where # stands for an integer either in the range
//                   [1,2147483647], or # = 0 which causes a seed to be
//                   chosen as the value of the UNIX clock in seconds; this
//                   guarantees that no two calls will give the same value for
//                   the seed if they are more than 2.0 seconds apart.
//option:   -m #
//                   mass fraction of the (infinitely extended) Plummer model
//option:   -r #
//                   radius fraction of the (infinitely extended) Plummer model
//
//            If the mass fraction and radius fraction are both 1.0 then 
//            particles will be sprinkled in all over space. If mfrac < 1.0 or
//            rfrac < 1.0 then each particle is constrained to lie within
//            both the radial and (cumulative) mass bound.
//               For example, if rfrac( mfrac ) > rfrac then rfrac is the 
//            limiting factor, but if rfrac( mfrac ) < rfrac then mfrac limits
//            the extent of the Plummer realization.
//
//option:   -i
//
//option:   -o
//
//option:   -c "..."
//
//-

/*
 *-----------------------------------------------------------------------------
 */
main(int argc, char ** argv)
    {
    int  i;
    int  c;
    int  n;
    int  input_seed, actual_seed;
    real  mfrac;
    real  rfrac;
    bool  c_flag = FALSE;
    bool  n_flag = FALSE;
    bool  s_flag = FALSE;
    bool  m_flag = FALSE;
    bool  r_flag = FALSE;
    bool  i_flag = FALSE;
    bool  o_flag = FALSE;
    char  *comment;
    char  seedlog[SEED_STRING_LENGTH];
    extern char *poptarg;
    int  pgetopt(int, char **, char *);

    if (argc == 1)
	{
	cerr <<"usage: mkplummer -n # [-s #] [-r #] [-m #] [-i] [-c \"..\"]\n";
	exit(1);
	}

    mfrac = rfrac = VERY_LARGE_NUMBER;

    while ((c = pgetopt(argc, argv, "c:n:s:m:r:io")) != -1)
	switch(c)
	    {
	    case 'n': n_flag = TRUE;
		      n = atoi(poptarg);
		      break;
	    case 's': s_flag = TRUE;
		      input_seed = atoi(poptarg);
		      break;
	    case 'm': m_flag = TRUE;
		      mfrac = atof(poptarg);
		      break;
	    case 'r': r_flag = TRUE;
		      rfrac = atof(poptarg);
		      break;
	    case 'i': i_flag = TRUE;
		      break;
	    case 'o': o_flag = TRUE;
		      break;
	    case 'c': c_flag = TRUE;
		      comment = poptarg;
		      break;
            case '?': cerr << "usage: mkplummer " <<
			      "-n # [-s #] [-r #] [-m #] [-i] [-c \"..\"]\n";
		      exit(1);
	    }            
    
    if (n_flag == FALSE)
	{
        cerr << "mkplummer: please specify the number # of";
	cerr << " particles with -n#\n";
	exit(1);
	}
//+todo:
// this should be moved elsewhere, of course.  Included here just to
// get it to work quickly (Piet, 921210).
//-todo:

    if (n < 1)
	{
        cerr << "mkplummer: n < 1 not allowed\n";
	exit(1);
        }

    dyn *b, *by, *bo;
    b = new dyn();
    bo = new dyn();
    if (i_flag)
        bo->set_label(1);
    b->set_oldest_daughter(bo);
    bo->set_parent(b);

    for (i = 1; i < n; i++)
        {
        by = new dyn();
	if (i_flag)
	    by->set_label(i+1);
	by->set_parent(b);
        bo->set_younger_sister(by);
        by->set_elder_sister(bo);
        bo = by;
        }

    if (c_flag == TRUE)
        b->log_comment(comment);
    b->log_history(argc, argv);

    if (s_flag == FALSE)
        input_seed = 0;                         /* default */
    actual_seed = srandinter(input_seed);

    if (o_flag) cerr << "mkplummer: random seed = " << actual_seed << endl;

    sprintf(seedlog, "         random number generator seed = %d",actual_seed);
    b->log_comment(seedlog);

    if (m_flag == FALSE && r_flag == FALSE)
        mfrac = MFRAC_DEFAULT;                    /* default */

    if (n != 0)
        mkplummer(b, n, mfrac, rfrac);

//  setup_starlab();
    cout.precision(2);
    put_dyn(cout, *b);
    }

#endif

/* end of: mkplummer.C */