#!/usr/bin/env python
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

nx, ny = (5,5)
Rf = 12
Rs = 12

x = np.linspace(0,5, nx*5)
y = np.linspace(0,5, ny*5)

xx, yy = np.meshgrid(x, y)
z = np.sqrt(Rf**2 - xx**2) - Rf + np.sqrt(Rs**2-yy**2) - Rs
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

h = ax.scatter(xx,yy,z)

f = open('testfile.asc', 'w')
xnew = np.dstack(xx)
xnew = np.reshape(xnew, -1)
ynew = np.dstack(yy)
znew = np.reshape(z,-1)
iter = 0
for inx in x:
	for iny in y:
		for inz in znew:
			#f.write('%d %d %d' % (inx, iny, inz))
			f.write('%f %f %f \n' % (inx, iny, znew[iter]))
			iter +=1
			break
#zneww = np.vstack(z)

print np.shape(znew)
print znew.size
print np.shape(x)

f.close()


#plt.show(h)