// generate n random particles orbiting a central one
// writes a particle2 file to output
// can be piped into nbody program, e.g. "orbit 1000 3 4 .2 .4 | nbody ."

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

static double G = 1;	// gravitational constant

double frand() {
    // return a random number between 0 and 1
    return rand()/32767.;
}

void main(int argc, char **argv) {
    if (argc!=6) {
	fprintf(stderr,
	    "Usage: orbit <sunmass> <#planets> <planetmass> <rmin> <rmax>\n");
	exit(1);
    }
    double sunmass = atof(argv[1]);
    int n = atoi(argv[2]);
    double planetmass = atof(argv[3]);
    double rmin = atof(argv[4]);
    double rmax = atof(argv[5]);
    int i;
    printf("particle2\n# orbit %g %d %g %g %g\n",
	sunmass, n, planetmass, rmin, rmax);
    printf("m %g\np .5 .5 0 0\nm %g\n", sunmass, planetmass);
    int seed = time(0);
    fprintf(stderr, "random seed=%d\n", seed);
    srand(seed);
    double speed = sqrt(G*sunmass);
	// assuming the potential is G*mi*mj*log(r)
    double r, theta;
    for (i=0; i<n; i++) {
	r = (rmax-rmin)*frand()+rmin;
	theta = 2*M_PI*frand();
	printf("p %g %g %g %g\n",
	    .5+r*cos(theta), .5+r*sin(theta),
	    -speed*sin(theta), speed*cos(theta));
    }
}