import sys
import re
import math
import time
import datetime
import numpy as np

import xml.etree.ElementTree as ET

f = open('globe.svg', 'w')
f.write('<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n<svg viewBox="0 0 800 800" version="1.1"\nxmlns="http://www.w3.org/2000/svg" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n')

tree = ET.parse('BlankMap-Equirectangular-states-3.svg')
root = tree.getroot()

camera = {'x': -15, 'y': 15, 'z': 30}
cameraForward = {'x': -1*camera['x'], 'y': -1*camera['y'], 'z': -1*camera['z']}
#cameraForward = {'x': -1, 'y': 3, 'z': 2}
cameraPerpendicular = {'x': cameraForward['y'], 'y': -1*cameraForward['x'], 'z': 0}
focalLength = 1000

radius = 10
pi = 3.14159

#magnitude of a 3D vector
def sumOfSquares(vector):
	return vector['x']**2 + vector['y']**2 + vector['z']**2
def magnitude(vector):
	return math.sqrt(sumOfSquares(vector))

#projects u onto v
def vectorProject(u, v):
	return np.dot(vectorToList (v), vectorToList (u))/magnitude(v)

#get unit vector
def unitVector(vector):
	magVector = magnitude(vector)
	return {'x': vector['x']/magVector, 'y': vector['y']/magVector, 'z': vector['z']/magVector }

#Calculates vector from two points, dictionary form
def findVector (origin, point):	
	return { 'x': point['x'] - origin['x'], 'y': point['y'] - origin['y'], 'z': point['z'] - origin['z'] }

#converts dictionary vector to list vector
def vectorToList (vector):
	return [vector['x'], vector['y'], vector['z']]

#Calculates horizon plane vector (points upward)
cameraHorizon = {'x': np.cross(vectorToList(cameraForward) , vectorToList(cameraPerpendicular))[0], 'y': np.cross(vectorToList(cameraForward) , vectorToList(cameraPerpendicular))[1], 'z': np.cross(vectorToList(cameraForward) , vectorToList(cameraPerpendicular))[2] }

#Takes a rectangular point and finds its height above the horizon (xDistance), displacement right or left of the perpendicular (yDistance) and z (zDistance). z is orthogonal distance to the plane facing the camera containing the point, not the point itself.
def physicalProjection (point):
	pointVector = findVector(camera, point)
		#pointVector is a vector starting from the camera and ending at a point in question
	return {'x': vectorProject(pointVector, cameraPerpendicular), 'y': vectorProject(pointVector, cameraHorizon), 'z': vectorProject(pointVector, cameraForward)}


# draws points onto camera sensor using xDistance, yDistance, and zDistance
def perspectiveProjection (pCoords):
	scaleFactor = focalLength/pCoords['z']
	return {'x': pCoords['x']*scaleFactor, 'y': pCoords['y']*scaleFactor}

#handles occlusion

cameraDistanceSquare = sumOfSquares(camera)
		#distance from globe center to camera

def distanceToPoint(spherePoint):
	point = sphereToRect(radius, spherePoint['long'], spherePoint['lat'])
	ray = findVector(camera,point)
	return vectorProject(ray, cameraForward)

def occlude(spherePoint):
	point = sphereToRect(radius, spherePoint['long'], spherePoint['lat'])
	ray = findVector(camera,point)
	d1 = magnitude(ray)
		#distance from camera to point

	dot_l = np.dot( [ray['x']/d1, ray['y']/d1, ray['z']/d1], vectorToList(camera) )
		#dot product of unit vector from camera to point and camera vector

	determinant = math.sqrt(abs( (dot_l)**2 - cameraDistanceSquare + radius**2 ))
	dNear = -(dot_l) + determinant
	dFar = -(dot_l) - determinant
		#formulas for finding the 2 points that the vector intersects the globe
	if d1 - 0.0000000001 <= dNear and d1 - 0.0000000001 <= dFar :
		return True
	else:
		return False

#converts spherical coordinates to rectangular coordinates
def sphereToRect (r, a, b):
	aRad = math.radians(a)
	bRad = math.radians(b)
	r_sin_b = r*math.sin(bRad)
	return {'x': r_sin_b*math.cos(aRad), 'y': r_sin_b*math.sin(aRad), 'z': r*math.cos(bRad) }
		#optimized
#	return {'x': r*math.sin(b*pi/180)*math.cos(a*pi/180), 'y': r*math.sin(b*pi/180)*math.sin(a*pi/180), 'z': r*math.cos(b*pi/180) }

#functions for plotting points
def rectPlot (coordinate):
	return perspectiveProjection(physicalProjection(coordinate))
def spherePlot (coordinate, sRadius):
	return rectPlot(sphereToRect(sRadius, coordinate['long'], coordinate['lat']))

#Draws SVG circle object
def svgCircle (coordinate, circleRadius, color):
	f.write('<circle cx=\"' + str(coordinate['x'] + 400) + '\" cy=\"' + str(coordinate['y'] + 400) + '\" r=\"' + str(circleRadius) + '\" style=\"fill:' + color + ';\"/>\n')
#Draws SVG polygon node
def polyNode (coordinate):
	f.write(str(coordinate['x'] + 400) + ',' + str(coordinate['y'] + 400) + ' ')

start_time = time.time()

#TRACE LIMB

rsquare_CDS = radius**2/cameraDistanceSquare
limbRadius = math.sqrt( radius**2 - radius**2*rsquare_CDS )

cx = camera['x']*rsquare_CDS
cy = camera['y']*rsquare_CDS
cz = camera['z']*rsquare_CDS

planeDistance = magnitude(camera)*(1 - rsquare_CDS)
planeDisplacement = math.sqrt(cx**2 + cy**2 + cz**2)

#trade & negate x and y to get a perpendicular vector
unitVectorCamera = unitVector(camera)
aV = unitVector( {'x': -unitVectorCamera['y'], 'y': unitVectorCamera['x'], 'z': 0} )
bV = np.cross(vectorToList(aV), vectorToList( unitVectorCamera ))

f.write('<polygon id=\"globe\" points=\"')
for t in range(180):
	theta = math.radians(2*t)
	cosT = math.cos(theta)
	sinT = math.sin(theta)
	
	limbPoint = { 'x': cx + limbRadius*(cosT*aV['x'] + sinT*bV[0]), 'y': cy + limbRadius*(cosT*aV['y'] + sinT*bV[1]), 'z': cz + limbRadius*(cosT*aV['z'] + sinT*bV[2]) }

	polyNode(rectPlot(limbPoint))

f.write('\" style="fill:#eee;stroke: none;stroke-width:0.5\" />')

f.write('<clipPath id=\"clipglobe\"><use xlink:href=\"#globe\"/></clipPath>')

#DRAW GRID

for x in range(72):
	for y in range(36):
		if occlude( { 'long': 5*x, 'lat': 5*y } ):
			svgCircle (spherePlot( { 'long': 5*x, 'lat': 5*y }, radius ), 1, '#fff')

#reads the svg map coordinates

tangentscale = (radius + planeDisplacement)/(pi*0.5)

miscDots = []
countries = len(root[4][0]) - 1
drawnCountries = 0
for x in range(countries):
	for path in root[4][0][x + 1]:
		for k in re.split('Z M', path.attrib['d']):
	
			countryIsVisible = False
			country = []
			for i in re.split('Z|M|L', k):
	
				breakup = re.split(',', re.sub("[^-0123456789.,]", "", i))
#				print (breakup)
				if breakup[0] and breakup[1]:
					sphereCoordinates = {}
					sphereCoordinates['long'] = float(breakup[0])
					sphereCoordinates['lat'] = float(breakup[1]) - 90
	
					#DRAW COUNTRY
	
					if occlude(sphereCoordinates):
						country.append([sphereCoordinates, radius])
	
						countryIsVisible = True
	
					else:
						rr = 1 + abs(math.tan( (distanceToPoint(sphereCoordinates) - planeDistance)/tangentscale ))
						country.append([sphereCoordinates, radius*rr])
	
			if countryIsVisible:
				f.write('<polygon clip-path="url(#clipglobe)" points=\"')
#				f.write('<polygon points=\"')
				for i in country:
					polyNode(spherePlot(i[0], i[1]))
				f.write('\" style="fill:#ff3092;stroke: #fff;stroke-width:0.3\" />\n\n')
				drawnCountries = drawnCountries + 1



f.write('</svg>')

elapsed_time = time.time() - start_time

print ('Drew ' + str(drawnCountries) + ' fragments from ' + str(countries) + ' countries in ' + str(elapsed_time) + ' seconds')