threebody.c

三体

/* -*- linux-c -*- */
/* threebody.c

   Copyright (C) 2002-2003 John M. Fregeau

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <getopt.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include "fewbody.h"
#include "threebody.h"

/* print the usage */
void print_usage(FILE *stream)
{
	fprintf(stream, "USAGE:\n");
	fprintf(stream, "  threebody [options...]\n");
	fprintf(stream, "\n");
	fprintf(stream, "OPTIONS:\n");
	fprintf(stream, "  -m --m0 <m0/MSUN>            : set mass of single star [%.6g]\n", M0/FB_CONST_MSUN);
	fprintf(stream, "  -n --m10 <m10/MSUN>          : set mass 0 of binary [%.6g]\n", M10/FB_CONST_MSUN);
	fprintf(stream, "  -o --m11 <m11/MSUN>          : set mass 1 of binary [%.6g]\n", M11/FB_CONST_MSUN);
	fprintf(stream, "  -a --a1 <a1/AU>                : set semimajor axis of binary [%.6g]\n", A1/FB_CONST_AU);
	fprintf(stream, "  -e --e1 <e1>                   : set eccentricity of binary 0 [%.6g]\n", E1);
	fprintf(stream, "  -v --vinf <vinf/v_crit>      : set velocity at infinity [%.6g]\n", VINF);
	fprintf(stream, "  -b --b <b/(a_0+a_1)>         : set impact parameter [%.6g]\n", B);
	fprintf(stream, "  -t --tstop <tstop/t_dyn>     : set stopping time [%.6g]\n", TSTOP);
	fprintf(stream, "  -D --dt <dt/t_dyn>           : set approximate output dt [%.6g]\n", DT);
	fprintf(stream, "  -c --tcpustop <tcpustop/sec> : set cpu stopping time [%.6g]\n", TCPUSTOP);
	fprintf(stream, "  -A --absacc <absacc>         : set integrator's absolute accuracy [%.6g]\n", ABSACC);
	fprintf(stream, "  -R --relacc <relacc>         : set integrator's relative accuracy [%.6g]\n", RELACC);
	fprintf(stream, "  -N --ncount <ncount>         : set number of integration steps between calls\n");
	fprintf(stream, "                                 to fb_classify() [%d]\n", NCOUNT);
	fprintf(stream, "  -z --tidaltol <tidaltol>     : set tidal tolerance [%.6g]\n", TIDALTOL);
	fprintf(stream, "  -k --ks                      : turn K-S regularization on or off [%d]\n", KS);
	fprintf(stream, "  -s --seed                    : set random seed [%ld]\n", SEED);
	fprintf(stream, "  -d --debug                   : turn on debugging\n");
	fprintf(stream, "  -V --version                 : print version info\n");
	fprintf(stream, "  -h --help                    : display this help text\n");
}

/* set up the scattering */
int init_scattering(fb_obj_t *obj[2], double vinf, double b, double tidaltol)
{
	double m0, m1, M, mu, a1, e1, m10, m11, r, x, y, vx, vy, qa, qb, qc;

	/* some useful variables */
	m0 = obj[0]->m;
	m1 = obj[1]->m;

	M = m0 + m1;
	mu = m0 * m1 / M;

	a1 = obj[1]->a;
	e1 = obj[1]->e;
	m10 = obj[1]->obj[0]->m;
	m11 = obj[1]->obj[1]->m;

	/* the separation that gives the desired tidal perturbation */
	r = pow(2.0*m0*m1/tidaltol, 1.0/3.0) * pow(m10*m11, -1.0/3.0)*a1*(1.0+e1);

	/* solve for y */
	qa = M*M + fb_sqr(b*vinf*vinf);
	qb = 2.0*(M*r-fb_sqr(b*vinf))*b*vinf*vinf;
	qc = -2.0*fb_sqr(b*vinf)*M*r + fb_sqr(fb_sqr(b*vinf));
	y = (-qb+sqrt(qb*qb-4.0*qa*qc))/(2.0*qa);

	/* x should never be zero, so no need to worry about the square root here */
	x = sqrt(r*r - y*y);

	/* determine v_x from quadratic */
	qa = 1.0 + fb_sqr(y/x);
	qb = 2.0*y*b*vinf/(x*x);
	qc = fb_sqr(b*vinf/x) - 2.0*M/r - vinf*vinf;
	vx = (-qb-sqrt(qb*qb-4.0*qa*qc))/(2.0*qa);

	vy = (b*vinf+y*vx)/x;

	/* set the positions and velocities */
	obj[0]->x[0] = - x / (1.0 + m0/m1);
	obj[0]->x[1] = - y / (1.0 + m0/m1);
	obj[0]->x[2] = 0.0;

	obj[0]->v[0] = - vx / (1.0 + m0/m1);
	obj[0]->v[1] = - vy / (1.0 + m0/m1);
	obj[0]->v[2] = 0.0;

	obj[1]->x[0] = x / (1.0 + m1/m0);
	obj[1]->x[1] = y / (1.0 + m1/m0);
	obj[1]->x[2] = 0.0;

	obj[1]->v[0] = vx / (1.0 + m1/m0);
	obj[1]->v[1] = vy / (1.0 + m1/m0);
	obj[1]->v[2] = 0.0;

	return(0);
}

/* calculate the units used */
int calc_units(fb_obj_t *obj[2], fb_units_t *units)
{
	units->v = sqrt(FB_CONST_G*(obj[0]->m + obj[1]->m)/(obj[0]->m * obj[1]->m) * \
			(obj[1]->obj[0]->m * obj[1]->obj[1]->m / obj[1]->a));
	units->l = obj[1]->a;
	units->t = units->l / units->v;
	units->m = units->l * fb_sqr(units->v) / FB_CONST_G;
	units->E = units->m * fb_sqr(units->v);

	return(0);
}

/* normalize all variables to our system of units */
int normalize(fb_obj_t *obj[2], fb_units_t units)
{
	int i, j, k;

	for (j=0; j<2; j++) {
		obj[j]->a /= units.l;
		obj[j]->m /= units.m;

		for (i=0; i<3; i++) {
			obj[j]->x[i] /= units.l;
			obj[j]->v[i] /= units.v;
		}
	}

	for (k=0; k<2; k++) {
		obj[1]->obj[k]->m /= units.m;

		for (i=0; i<3; i++) {
			obj[1]->obj[k]->x[i] /= units.l;
			obj[1]->obj[k]->v[i] /= units.v;
		}
	}

	return(0);
}

/* the main attraction */
int main(int argc, char *argv[])
{
	int i, j;
	unsigned long int seed;
	double m0, m10, m11, a1, e1;
	double vinf, b, m1, M, mu, Ei, E, Li[3], l0[3], l1[3], L[3], r[3], t;
	fb_hier_t hier;
	fb_input_t input;
	fb_ret_t retval;
	fb_units_t units;
	char string1[1024], string2[1024];
	gsl_rng *rng;
	const gsl_rng_type *rng_type=gsl_rng_mt19937;
	const char *short_opts = "m:n:o:a:e:v:b:t:D:c:A:R:N:z:k:s:dVh";
	const struct option long_opts[] = {
		{"m0", required_argument, NULL, 'm'},
		{"m10", required_argument, NULL, 'n'},
		{"m11", required_argument, NULL, 'o'},
		{"a1", required_argument, NULL, 'a'},
		{"e1", required_argument, NULL, 'e'},
		{"vinf", required_argument, NULL, 'v'},
		{"b", required_argument, NULL, 'b'},
		{"tstop", required_argument, NULL, 't'},
		{"dt", required_argument, NULL, 'D'},
		{"tcpustop", required_argument, NULL, 'c'},
		{"absacc", required_argument, NULL, 'A'},
		{"relacc", required_argument, NULL, 'R'},
		{"ncount", required_argument, NULL, 'N'},
		{"tidaltol", required_argument, NULL, 'z'},
		{"ks", required_argument, NULL, 'k'},
		{"seed", required_argument, NULL, 's'},
		{"debug", no_argument, NULL, 'd'},
		{"version", no_argument, NULL, 'V'},
		{"help", no_argument, NULL, 'h'},
		{NULL, 0, NULL, 0}
	};

	/* set parameters to default values */
	m0 = M0;
	m10 = M10;
	m11 = M11;
	a1 = A1;
	e1 = E1;
	vinf = VINF;
	b = B;
	input.ks = KS;
	input.tstop = TSTOP;
	input.Dflag = 0;
	input.dt = DT;
	input.tcpustop = TCPUSTOP;
	input.absacc = ABSACC;
	input.relacc = RELACC;
	input.ncount = NCOUNT;
	input.tidaltol = TIDALTOL;
	seed = SEED;
	fb_debug = DEBUG;

	while ((i = getopt_long(argc, argv, short_opts, long_opts, NULL)) != -1) {
		switch (i) {
		case 'm':
			m0 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'n':
			m10 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'o':
			m11 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'a':
			a1 = atof(optarg) * FB_CONST_AU;
			break;
		case 'e':
			e1 = atof(optarg);
			if (e1 >= 1.0) {
				fprintf(stderr, "e0 must be less than 1\n");
				return(1);
			}
			break;
		case 'v':
			vinf = atof(optarg);
			if (vinf < 0.0) {
				fprintf(stderr, "vinf must be non-negative\n");
				return(1);
			}
			break;
		case 'b':
			b = atof(optarg);
			if (b < 0.0) {
				fprintf(stderr, "b must be non-negative\n");
				return(1);
			}
			break;
		case 't':
			input.tstop = atof(optarg);
			break;
		case 'D':
			input.Dflag = 1;
			input.dt = atof(optarg);
			break;
		case 'c':
			input.tcpustop = atof(optarg);
			break;
		case 'A':
			input.absacc = atof(optarg);
			break;
		case 'R':
			input.relacc = atof(optarg);
			break;
		case 'N':
			input.ncount = atoi(optarg);
			break;
		case 'z':
			input.tidaltol = atof(optarg);
			break;
		case 'k':
			input.ks = atoi(optarg);
			break;
		case 's':
			seed = atol(optarg);
			break;
		case 'd':
			fb_debug = 1;
			break;
		case 'V':
			fb_print_version(stdout);
			return(0);
		case 'h':
			fb_print_version(stdout);
			fprintf(stdout, "\n");
			print_usage(stdout);
			return(0);
		default:
			break;
		}
	}

	/* check to make sure there was nothing crazy on the command line */
	if (optind < argc) {
		print_usage(stdout);
		return(1);
	}

	/* initialize a few things for integrator */
	t = 0.0;
	hier.nstar = 3;
	fb_malloc_hier(&hier);
	fb_init_hier(&hier);

	/* put stuff in log entry */
	sprintf(&(input.firstlogentry[strlen(input.firstlogentry)]), "  command line:");
	for (i=0; i<argc; i++) {
		sprintf(&(input.firstlogentry[strlen(input.firstlogentry)]), " %s", argv[i]);
	}
	sprintf(&(input.firstlogentry[strlen(input.firstlogentry)]), "\n");

	/* print out values of paramaters */
	fprintf(stderr, "PARAMETERS:\n");
	fprintf(stderr, "  ks=%d  seed=%ld\n", input.ks, seed);
	fprintf(stderr, "  m0=%.6g MSUN\n", m0/FB_CONST_MSUN);
	fprintf(stderr, "  a1=%.6g AU  e1=%.6g  m10=%.6g MSUN  m11=%.6g MSUN\n", \
		a1/FB_CONST_AU, e1, m10/FB_CONST_MSUN, m11/FB_CONST_MSUN);
	fprintf(stderr, "  vinf=%.6g  b=%.6g  tstop=%.6g  tcpustop=%.6g\n", \
		vinf, b, input.tstop, input.tcpustop);
	fprintf(stderr, "  tidaltol=%.6g  abs_acc=%.6g  rel_acc=%.6g  ncount=%d\n\n", \
		input.tidaltol, input.absacc, input.relacc, input.ncount);

	/* initialize GSL rng */
	gsl_rng_env_setup();
	rng = gsl_rng_alloc(rng_type);
	gsl_rng_set(rng, seed);

	/* create binary */
	hier.hier[hier.hi[2]+0].obj[0] = &(hier.hier[hier.hi[1]+1]);
	hier.hier[hier.hi[2]+0].obj[1] = &(hier.hier[hier.hi[1]+2]);
	hier.hier[hier.hi[2]+0].t = t;

	/* give the objects some properties */
	for (j=0; j<hier.nstar; j++) {
		hier.hier[hier.hi[1]+j].id = j;
		sprintf(hier.hier[hier.hi[1]+j].idstring, "%d", j);
		hier.hier[hier.hi[1]+j].n = 1;
		hier.hier[hier.hi[1]+j].obj[0] = NULL;
		hier.hier[hier.hi[1]+j].obj[1] = NULL;
	}

	hier.hier[hier.hi[1]+0].m = m0;
	hier.hier[hier.hi[1]+1].m = m10;
	hier.hier[hier.hi[1]+2].m = m11;

	hier.hier[hier.hi[2]+0].m = m10 + m11;

	hier.hier[hier.hi[2]+0].a = a1;
	hier.hier[hier.hi[2]+0].e = e1;

	hier.obj[0] = &(hier.hier[hier.hi[1]+0]);
	hier.obj[1] = &(hier.hier[hier.hi[2]+0]);
	hier.obj[2] = NULL;

	/* get the units and normalize */
	calc_units(hier.obj, &units);
	normalize(hier.obj, units);

	fprintf(stderr, "UNITS:\n");
	fprintf(stderr, "  v=v_crit=%.6g km/s  l=%.6g AU  t=t_dyn=%.6g yr\n", \
		units.v/1.0e5, units.l/FB_CONST_AU, units.t/FB_CONST_YR);
	fprintf(stderr, "  M=%.6g M_sun  E=%.6g erg\n\n", units.m/FB_CONST_MSUN, units.E);

	fb_dprintf("moving hierarchies analytically in from infinity...\n");

	m0 = hier.obj[0]->m;
	m1 = hier.obj[1]->m;
	M = m0 + m1;
	mu = m0 * m1 / M;

	Ei = 0.5 * mu * fb_sqr(vinf);

	/* move hierarchies analytically in from infinity along hyperbolic orbit */
	init_scattering(hier.obj, vinf, b, input.tidaltol);

	/* and check to see that we conserved energy and angular momentum */
	fb_cross(hier.obj[0]->x, hier.obj[0]->v, l0);
	fb_cross(hier.obj[1]->x, hier.obj[1]->v, l1);

	for (i=0; i<3; i++) {
		L[i] = (m0 * l0[i] + m1 * l1[i]);
		r[i] = hier.obj[1]->x[i] - hier.obj[0]->x[i];
	}

	E = - m0 * m1 / fb_mod(r) + 0.5 * (m0 * fb_dot(hier.obj[0]->v, hier.obj[0]->v) + \
					m1 * fb_dot(hier.obj[1]->v, hier.obj[1]->v));

	fb_dprintf("L0=%.6g DeltaL/L0=%.6g DeltaL=%.6g\n", mu*b*vinf, fb_mod(L)/(mu*b*vinf)-1.0, fb_mod(L)-mu*b*vinf);
	fb_dprintf("E0=%.6g DeltaE/E0=%.6g DeltaE=%.6g\n\n", Ei, E/Ei-1.0, E-Ei);

	/* trickle down the binary properties, then back up */
	fb_dprintf("obj[%d]->a=%e\n", 1, hier.obj[1]->a);
	fb_dprintf("obj[%d]->e=%e\n", 1, hier.obj[1]->e);
	fb_dprintf("obj[%d]->m=%e\n", 1, hier.obj[1]->m);
	fb_randorient(&(hier.hier[hier.hi[2]+0]), rng);
	fb_downsync(&(hier.hier[hier.hi[2]+0]), t);
	fb_upsync(&(hier.hier[hier.hi[2]+0]), t);
	fb_dprintf("obj[%d]->a=%e\n", 1, hier.obj[1]->a);
	fb_dprintf("obj[%d]->e=%e\n", 1, hier.obj[1]->e);
	fb_dprintf("obj[%d]->m=%e\n", 1, hier.obj[1]->m);

	/* store the initial energy and angular momentum*/
	Ei = fb_petot(&(hier.hier[hier.hi[1]]), hier.nstar) + fb_ketot(&(hier.hier[hier.hi[1]]), hier.nstar);
	fb_angmom(&(hier.hier[hier.hi[1]]), hier.nstar, Li);

	/* integrate along */
	fb_dprintf("calling fewbody()...\n");

	/* call fewbody! */
	retval = fewbody(input, &hier, &t);

	/* print information to screen */
	fprintf(stderr, "OUTCOME:\n");
	if (retval.retval == 1) {
		fprintf(stderr, "  encounter complete: %s (%s)\n\n", \
			fb_sprint_hier(hier.obj, hier.nobj, string1), \
			fb_sprint_hier_hr(hier.obj, hier.nstar, hier.nobj, string2));
	} else {
		fprintf(stderr, "  encounter NOT complete: %s (%s)\n\n", \
			fb_sprint_hier(hier.obj, hier.nobj, string1), \
			fb_sprint_hier_hr(hier.obj, hier.nstar, hier.nobj, string2));
	}

	fb_dprintf("there were %ld integration steps\n", retval.count);
	fb_dprintf("fb_classify() was called %ld times\n", retval.iclassify);

	fprintf(stderr, "FINAL:\n");
	fprintf(stderr, "  t_final=%.6g (%.6g yr)  t_cpu=%.6g s\n", \
		t, t*units.t/FB_CONST_YR, retval.tcpu);

	fprintf(stderr, "  L0=%.6g  DeltaL/L0=%.6g  DeltaL=%.6g\n", fb_mod(Li), retval.DeltaLfrac, retval.DeltaL);
	fprintf(stderr, "  E0=%.6g  DeltaE/E0=%.6g  DeltaE=%.6g\n", Ei, retval.DeltaEfrac, retval.DeltaE);
	fprintf(stderr, "  Rmin=%.6g (%.6g RSUN)  Rmin_i=%d  Rmin_j=%d\n", \
		retval.Rmin, retval.Rmin*units.l/FB_CONST_RSUN, retval.Rmin_i, retval.Rmin_j);

	/* free GSL stuff */
	gsl_rng_free(rng);

	/* free our own stuff */
	fb_free_hier(hier);

	/* done! */
	return(0);
}