#include <iostream>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include "movie.h"
#include "putable.h"
using namespace std;

/*
Radius of Milky Way is 50,000 light years.
It will be centered at the origin and lie in the x-y plane.
Its radius will be picture::ymax - 5.
Length of bar is 27,000 light years.
Its radius will be (picture::ymax - 5) * 13,500 / 50,000.

Four arms of Milky Way are logarithmic spirals:
	r = a * pow(b, theta)
	  = a * exp(theta * lnb).
theta = log to the base b of (r/a) = (ln r - ln a) / ln b.

Pitch of arms is 12 degrees.
ln b = tan (12 * pi / 180), so b = exp(tan(12 * pi / 180)).
Arms extend from end of bar, so a = (picture::ymax - 5) * 13,000 / 50,000.
Therefore theta will range from 0 to (ln (picture::ymax - 5) - ln a) / ln b.
*/

int main(int argc, char **argv)
{
	srand(time(0));

	const double lnb = tan(12 * pi / 180);	//natural log of b
	const double b = exp(lnb);
	const double radius = picture::ymax - 5; 	//of galaxy
	const double a = radius * 13500 / 50000;
	const double theta_max = (log(radius) - log(a)) / log(b);
	
	const size_t narms = 4;			//number of arms
	const size_t npoints = 220;		//number of points per arm
	const size_t npbar = 9 * npoints / 10;	//number of points in bar
	const size_t n = npbar + npoints * narms;

	//two-thirds the interval between arms at the armstart
	const double interval = 2 * a * (pow(b, pi / 2) - 1) / 3;
	xyz arr[n];

	for (size_t i = 0; i < npbar; ++i) {
		double x = a * myrand();
		if (rand() % 2) {
			x = -x;
		}

		double y = interval * weighted_rand();
		if (rand() % 2) {
			y = -y;
		}

		double z = a * a * a / (2 * (x * x + a * a));
		z *= weighted_rand();
		if (rand() % 2) {
			z = -z;
		}
		arr[i] = xyz(x, y, z);
	}

	for (size_t i = 0; i < npoints; ++i) {
		for (size_t arm = 0; arm < narms; ++arm) {
			const double theta = theta_max * myrand();
			double r = a * exp(theta * lnb);
			double dr = interval;
			dr *= weighted_rand();
			if (rand() % 2) {
				dr = -dr;
			}
			r += dr;

			//Witch of Agnesi
			double dz = a * a * a / (2 * (r * r + a * a));
			dz *= weighted_rand();
			if (rand() % 2) {
				dz = -dz;
			}
			xyz point(r * cos(theta), r * sin(theta), dz);

			point.zrot(arm * pi / 2);
			arr[npbar + i + arm * npoints] = point;
			/*
			cerr << "arr[" << npbar + i + arm * npoints << "] == "
				<< point << "\n";
			*/
		}
	}
	putable stars;
	stars.putable_set(arr, arr + n);

	const picture empty;

	const xyz viewpoint(0, -picture::dmax);
	picture fade_in_milky(viewpoint, xyz());
	fade_in_milky << stars;

	movie(empty, fade_in_milky).write(
		"fade_in_milky.gif",
		"empty",
		"Milky Way edge-on");

	movie milky_to_overhead(
		xarc(viewpoint, xyz(0, 0, picture::dmax)),
		constant(xyz()),
		constant(picture::dmax),
		3 * movie::nf / 2
	);
	(milky_to_overhead << stars).write(
		"milky_to_overhead.gif",
		"Milky Way edge-on",
		"Milky Way from above"
	);

	system(
		"rm -f milky_to_edge.gif;"
		"ln milky_to_overhead_r.gif milky_to_edge.gif;"

		"rm -f milky_to_edge_r.gif;"
		"ln milky_to_overhead.gif milky_to_edge_r.gif;"

		"rm -f fade_out_milky.gif;"
		"ln fade_in_milky_r.gif fade_out_milky.gif;"

		"rm -f fade_out_milky_r.gif;"
		"ln fade_in_milky.gif fade_out_milky_r.gif;"
	);

	return EXIT_SUCCESS;
}