/* sphere - creates a svg vector-graphics file which depicts a wireframe sphere
 *
 * Copyright (C) 2008 Wikimedia foundation
 *
 * 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, 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, you can either send email to this
 * program's author (see below) or write to:
 *   The Free Software Foundation, Inc.
 *   51 Franklin Street, Fifth Floor
 *   Boston, MA 02110-1301  USA
 */

/* The expressions in this code are not proven to be correct.
 * Hence this code probably contains lots of bugs. Be aware! */

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>

using namespace std;

const double PI = 3.1415926535897932;
const double DEG = PI / 180.0;

/********************************* settings **********************************/
int n_lon = 18; 			// number of latitude fields (18 => 10° each)
int n_lat = 18; 			// half number of longitude fields (18 => 10° each)
double lon_offset = 2.5 * DEG; 	// offset of the meridians
double w = 52.5 * DEG; 		// axial tilt (0° => axis is perpendicular to image plane)
double stripe_grad = 0.5 * DEG;	// width of each line
int image_size = 400;			// width and height of the image in pixels
double back_opacity = 0.25;		// opacity of the sphere's backside
char color[] = "#334070";		// color of lines
int istep = 2; 			// svg code indentation step
/*****************************************************************************/

double sqr(double x)
{
	return(x * x);
}

// commands for svg-code:
void indent(int n, bool in_tag = false)
{
	n *= istep;
	if (in_tag) n += istep + 1;
	for (int i = 0; i < n; i++) cout << " ";
}
void M()
{
	cout << "M ";
}
void Z()
{
	cout << "Z ";
}
void xy(double x, double y)
{
	cout << x << ",";
	cout << y << " ";
}
void arc(double a, double b, double x_axis_rot, bool large_arc, bool sweep)
{	// draws an elliptic arc
	if (b < 0.5E-6)
	{	// flat ellipses are not rendered properly => use line
		cout << "L ";
	}
	else
	{
		cout << "A ";
		cout << a << ",";	// semi-major axis
		cout << b << " ";	// semi-minor axis
		cout << x_axis_rot << " ";
		cout << large_arc << " ";
		cout << sweep << " ";
	}
}
void circle(bool clockwise)
{
	M();
	xy(-1, 0);
	arc(1, 1, 0, 0, !clockwise);
	xy(1, 0);
	arc(1, 1, 0, 0, !clockwise);
	xy(-1, 0);
	Z();
}

void start_svg_file()
{
	cout << "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n";
	cout << "<svg id=\"Sphere_wireframe\"\n";
	cout << "  version=\"1.1\"\n";
	cout << "  baseProfile=\"full\"\n";
	cout << "  xmlns=\"http://www.w3.org/2000/svg\"\n";
	cout << "  xmlns:xlink=\"http://www.w3.org/1999/xlink\"\n";
	cout << "  width=\"" << image_size << "\"\n";
	cout << "  height=\"" << image_size << "\">\n\n";
	cout << "  <title>Sphere wireframe</title>\n\n";
	cout << "  <desc>\n";
	cout << "     about: http://commons.wikimedia.org/wiki/Image:Sphere_wireframe.svg\n";
	cout << "     rights: GNU Free Documentation license,\n";
	cout << "             Creative Commons Attribution ShareAlike license\n";
	cout << "  </desc>\n\n";
	cout << "  <g id=\"sphere\" transform=\"scale(" << 0.5 * image_size;
	cout << ", " << -0.5 * image_size << ") translate(1, -1)\">\n";
}

void end_svg_file()
{
	cout << "  </g>\n</svg>\n";
}



int main (int argc, char *argv[])
{
	// accept -lat and -lon as parameter
	for (int i = 2; i < argc; i++)
	{
		if (isdigit(argv[i][0]) || (sizeof(argv[i]) > sizeof(char)
			&& isdigit(argv[i][1])
			&& (argv[i][0] == '.' || argv[i][0] == '-')))
		{
			if (strcmp(argv[i - 1], "-lon") == 0)
			{
				lon_offset = atof(argv[i]) * DEG;
			}
			if (strcmp(argv[i - 1], "-lat") == 0)
			{
				w = atof(argv[i]) * DEG;
			}
		}
	}
	double cosw = cos(w), sinw = sin(w);
	double d = 0.5 * stripe_grad;

	start_svg_file();
	int ind = 2; // initial indentation level
	indent(ind);
	cout << "<g id=\"sphere_back\" transform=\"rotate(180)\" ";
	cout << "opacity=\"" << back_opacity << "\">\n";
	indent(++ind);
	cout << "<g id=\"sphere_half\">\n";

	// meridians
	indent(++ind); cout << "<g id=\"meridians\"\n";
	indent(ind++, true);
	cout << "style=\"stroke:none; fill:" << color << "; fill_rule:evenodd\">\n";
	double a = abs(cos(d));
	for (int i_lon = 0; i_lon < n_lat; i_lon++)
	{	// draw one meridian
		double longitude = lon_offset + (i_lon * 180.0 / n_lat) * DEG;
		double lon[2];
		lon[0] = longitude + d;
		lon[1] = longitude - d;
		
		indent(ind);
		cout << "<path id=\"meridian";
		cout << i_lon << "\"\n";
		indent(ind, true);
		cout << "d=\"";

		double axis_rot = atan2(-1.0 / tan(longitude), cosw);
		if (sinw < 0)
			axis_rot += PI;
		double w2 = sin(longitude) * sinw;
		double b = abs(w2 * cos(d));

		double sinw1 = sin(d) / sqrt(1.0 - sqr(sin(longitude) * sinw));

		if (abs(sinw1) >= 1.0)
		{	// stripe covers edge of the circle
			double w3 = sqrt(1.0 - sqr(w2)) * sin(d);
			circle(false);
			// ellipse
			M();
			xy(sin(axis_rot) * w3 - cos(axis_rot) * a,
				-cos(axis_rot) * w3 - sin(axis_rot) * a);
			arc(a, b, axis_rot / DEG, 0, 0);
			xy(sin(axis_rot) * w3 + cos(axis_rot) * a,
				-cos(axis_rot) * w3 + sin(axis_rot) * a);
			arc(a, b, axis_rot / DEG, 0, 0);
			xy(sin(axis_rot) * w3 - cos(axis_rot) * a,
				-cos(axis_rot) * w3 - sin(axis_rot) * a);
			Z();
		}
		else
		{	// draw a disrupted ellipse bow
			double w1 = asin(sinw1);
			M();
			xy(-cos(axis_rot + w1), -sin(axis_rot + w1));
			arc(a, b, axis_rot / DEG, 1, 0);
			xy(cos(axis_rot - w1), sin(axis_rot - w1));
			arc(1, 1, 0, 0, 1);
			xy(cos(axis_rot + w1), sin(axis_rot + w1));
			arc(a, b, axis_rot / DEG, 0, 1);
			xy(-cos(axis_rot - w1), -sin(axis_rot - w1));
			arc(1, 1, 0, 0, 1);
			xy(-cos(axis_rot + w1), -sin(axis_rot + w1));
		}
		Z();
		cout << "\" />\n";
	}
	indent(--ind); cout << "</g>\n";

	cout << endl;

	// circles of latitude
	indent(ind); cout << "<g id=\"circles_of_latitude\"\n";
	indent(ind, true);
	cout << "style=\"stroke:none; fill:" << color << "; fill_rule:evenodd\">\n";
	ind++;
	for (int i_lat = 1; i_lat < n_lon; i_lat++)
	{	// draw one circle of latitude
		double latitude = (i_lat * 180.0 / n_lon - 90.0) * DEG;
		double lat[2];
		lat[0] = latitude + d;
		lat[1] = latitude - d;
		double x[2], yd[2], ym[2];
		for (int i = 0; i < 2; i++)
		{
			x[i] = abs(cos(lat[i]));
			yd[i] = abs(cosw * cos(lat[i]));
			ym[i] = sinw * sin(lat[i]);
		}
		double h[4];	// height of each point above image plane
		h[0] = sin(lat[0] + w);
		h[1] = sin(lat[0] - w);
		h[2] = sin(lat[1] + w);
		h[3] = sin(lat[1] - w);
		
		if (h[0] > 0 || h[1] > 0 || h[2] > 0 || h[3] > 0)
		{	// at least any part visible
			indent(ind);
			cout << "<path id=\"circle_of_latitude";
			cout << i_lat << "\"\n";
			indent(ind, true);
			cout << "d=\"";
			for (int i = 0; i < 2; i++)
			{
				if ((h[2*i] >= 0 && h[2*i+1] >= 0)
					&& (h[2*i] > 0 || h[2*i+1] > 0))
				{	// complete ellipse
					M();
					xy(-x[i], ym[i]); // startpoint
					for (int z = 1; z > -2; z -= 2)
					{
						arc(x[i], yd[i], 0, 1, i);
						xy(z * x[i], ym[i]);
					}
					Z();
					if (h[2-2*i] * h[3-2*i] < 0)
					{	// partly ellipse + partly circle
						double yp = sin(lat[1-i]) / sinw;
						double xp = sqrt(1.0 - sqr(yp));
						if (sinw < 0)
						{
							xp = -xp;
						}
						M();
						xy(-xp, yp);
						arc(x[1-i], yd[1-i], 0,
							sin(lat[1-i]) * cosw > 0, cosw >= 0);
						xy(xp, yp);
						arc(1, 1, 0, 0, cosw >= 0);
						xy(-xp, yp);
						Z();
					}
					else if (h[2-2*i] <= 0 && h[3-2*i] <= 0)
					{	// stripe covers edge of the circle
						circle(cosw < 0);
					}
				}
			}
			
			if ((h[0] * h[1] < 0 && h[2] <= 0 && h[3] <= 0)
				|| (h[0] <= 0 && h[1] <= 0 && h[2] * h[3] < 0))
			{
				// one slice visible
				int i = h[0] <= 0 && h[1] <= 0;
				double yp = sin(lat[i]) / sinw;
				double xp = sqrt(1.0 - yp * yp);
				M();
				xy(-xp, yp);
				arc(x[i], yd[i], 0, sin(lat[i]) * cosw > 0, cosw * sinw >= 0);
				xy(xp, yp);
				arc(1, 1, 0, 0, cosw * sinw < 0);
				xy(-xp, yp);
				Z();
			}
			else if (h[0] * h[1] < 0 && h[2] * h[3] < 0)
			{
				// disrupted ellipse bow
				double xp[2], yp[2];
				for (int i = 0; i < 2; i++)
				{
					yp[i] = sin(lat[i]) / sinw;
					xp[i] = sqrt(1.0 - sqr(yp[i]));
					if (sinw < 0) xp[i] = -xp[i];
				}
				M();
				xy(-xp[0], yp[0]);
				arc(x[0], yd[0], 0, sin(lat[0]) * cosw > 0, cosw >= 0);
				xy(xp[0], yp[0]);
				arc(1, 1, 0, 0, 0);
				xy(xp[1], yp[1]);
				arc(x[1], yd[1], 0, sin(lat[1]) * cosw > 0, cosw < 0);
				xy(-xp[1], yp[1]);
				arc(1, 1, 0, 0, 0);
				xy(-xp[0], yp[0]);
				Z();
			}
			cout << "\" />\n";
		}
	}
	for (int i = 0; i < 3; i++)
	{
		indent(--ind);
		cout << "</g>\n";
	}
	indent(ind--);
	cout << "<use id=\"sphere_front\" xlink:href=\"#sphere_half\" />\n";
	end_svg_file();
}