import numpy as np
import matplotlib.pyplot as plt

# Parabolic cooker generator
# Written by John Wallin 2019


# the width of the parabola (inches)
xmax = 11.5

# the focal length (inches)
qq = 12

# number of points in the parabola
n = 100  

# allocate the lists needed for the points
xx =[]
yy = []
yy1 = []
plength = []
dy = []


# R is the focal point of the cooker - 2 times the focal length
R = qq * 2 

# A helper variable for the spatial step size
dx = xmax / float(n-1)


intY = 0.0
for i in range(n):

    # find the xposition
    x = float(i ) * dx

    # calculate y using the parabola formula 
    y = x*x / (4*qq)

    # save the results to a list
    xx.append(x)
    yy.append(y)

    # parabolic length - determined by doing a path integral along the curve
    dydx = x / (2 * qq)
    intY = intY + np.sqrt( 1.0 + dydx*dydx) * dx
    plength.append(intY)

    # calculate the spherical approximation of the curve
    y1 = R - np.sqrt( R*R- x*x)
    yy1.append(y1)
    
    # calcualte the difference between the sphere and the parabola 
    dy.append(y1-y)

    # the next line can be uncommented if you want to look at the data generated
    print x, y, y1, R*R, x*x, intY

# summarize the results"
print "focal length = ", qq
print "focal point = ", R
print "mirror length = ", intY *2
print "depth of mirror = ", yy[-1] - yy[0]
print "maximum deviation from a parabola = ", yy1[-1]- yy[-1]

# if you just want to look at the parabola, set showparabola = 1
# this will show the parabola and the circular approximation
showparabola = 1
if showparabola == 1:
    plt.plot(xxp2, yyp2, '-b')
    plt.plot(xxp1, yyp, '-b')
    plt.show()
    exit()

# set up the parameters for the svg plotting file
pwidth =  15 
pheight = 24

# initialize the plot routine
fig = plt.figure(figsize=(pwidth*2, pheight)) 
ax1 = fig.add_subplot(111)
ax1.axis("equal")
pwidth1 = 13
ax1.set_xlim([-pwidth,pwidth])
ax1.set_ylim([0,pheight])
plt.axis('off')

# set the parameters for the tabs on the sizes of the plots
yvoffset = 2.
yoffset = 1.
sidetabwidth = 1.5
sidetabbuffer = 0.5

# initialize helper arrays
xmaster = []
ymaster = []
yyp1= []
xxp1 = []
xxp2 = []
yyp2 = []

# since the parabola we calculated started at zero and went to xmax,
# we need to reflect that.  We do this such that the parabola and the 
# CNC path will be a simple closed path.
for i in range(len(xx)):

    # x goes from zero to xmax
    x = xx[i]
    y = yy[i]
    y = y + yoffset
    xxp1.append(x)
    yyp1.append(y)

    # x goes form -xmax to zero
    j = len(xx) - 1 - i
    x = xx[j]
    xxp2.append(-x)
    y = yy[j]
    y = y + yoffset
    yyp2.append(y)


# add all the parabolic points from -xmax to zero
xmaster = xmaster + xxp2
ymaster = ymaster + yyp2

# add the remaining parabolic points from zero to xmax
xmaster = xmaster + xxp1
ymaster = ymaster + yyp1


# set the left boundary
xleft = xxp2[0]
xright = xxp1[-1]

# set the tab boundary
xxleft = xxp2[0] - sidetabwidth 
xxright = xxp1[-1] + sidetabwidth


# this set of code systematically adds the line segments to make the tabs for the cuts
doconvex = 1
ytop =  yvoffset + yyp1[-1]
ybot = yyp1[0] - 0.5
if doconvex:

    # 2nd horizontal - right
    xmaster = xmaster + [xright, xright+sidetabbuffer]
    ymaster = ymaster + [yyp1[-1], yyp1[-1]]

    # 2nd vertical - right
    xmaster = xmaster + [xright+sidetabbuffer, xright+sidetabbuffer]
    ymaster = ymaster + [yyp1[-1],ybot]

    # bottom horizontal - right
    xmaster = xmaster + [xright+sidetabbuffer,xxright] 
    ymaster = ymaster + [ybot,ybot]

    # verticals- right
    xmaster = xmaster + [xxright,xxright]
    ymaster = ymaster + [ybot,ytop]

    # top horizontal
    xmaster = xmaster + [xxright,xxleft] 
    ymaster = ymaster + [ytop,ytop] 
    
    # vertical left
    xmaster = xmaster + [xxleft,xxleft] 
    ymaster = ymaster + [ytop,ybot] 

    # bottom horizontal - left
    xmaster = xmaster + [xxleft,xleft-sidetabbuffer] 
    ymaster = ymaster + [ybot,ybot] 

    # 2nd vertical - left
    #plt.plot( [xleft-sidetabbuffer,xleft-sidetabbuffer] , [ybot,yyp2[0]] , '-b')
    xmaster = xmaster + [xleft-sidetabbuffer,xleft-sidetabbuffer] 
    ymaster = ymaster + [ybot,yyp2[0]] 

    # 2nd horizontal - left
    #plt.plot( [xleft-sidetabbuffer,xleft] , [yyp2[0],yyp2[0]] , '-b')
    xmaster = xmaster + [xleft-sidetabbuffer,xleft] 
    ymaster = ymaster + [yyp2[0],yyp2[0]]  
        


        
# actually plot the path
plt.plot( xmaster, ymaster, '-r')        

plt.gca().set_position([0, 0, 1, 1])
plt.savefig("sundogger.svg")

# To import this file into easle, you need to do some minor edits to the SVG file.
# This is the diffs between the edits I did manually.   The 11a12,19 refers to lines 
# 12 to 19 in the original file being deleted.  

#11a12,19
#>   <g id="patch_1">
#>    <path d="M 0 1728 
#> L 2160 1728 
#> L 2160 0 
#> L 0 0 
#> z
#> " style="fill:#ffffff;"/>
#>   </g>
#14c22
#<     <path d="M 252 1457.625 
#---
#>     <path clip-path="url(#p32e7fa2eaf)" d="M 252 1457.625 
#112a121,122
#> L 1908 1457.625 
#> L 1944 1457.625 
#114a125,126
#> L 1944 1692 
#> L 2016 1692 
#116a129
#> L 2016 1313.625 
#117a131,132
#> L 144 1313.625 
#> L 144 1692 
#119a135
#> L 216 1692 
#120a137,138
#> L 216 1457.625 
#> L 252 1457.625 
#125a144,148
#>  <defs>
##>   <clipPath id="p32e7fa2eaf">
#>    <rect height="1728" width="2160" x="0" y="0"/>
#>   </clipPath>
#>  </defs>