reflectSurface = function(x, xline, fx, fxline, fpxline){



slope = fpxline
parab = fx
parabline = fxline

intercept = parabline-slope*xline

#reflections

mp = tan(2*atan(slope)-pi/2)
bp = parabline-mp*xline

plot(xline, parabline, t="p", lwd=2)
plot(x, parab, t="l", lwd=1)

for(xi in 1:length(xline)){
    #vertical rays
    abline(v=xline[xi])
    #tangents
    #abline(a=intercept[xi], b=slope[xi])
    #reflections
    abline(a=bp[xi], b=mp[xi])
    #points(xline[xi], parabline[xi], pch=23)
    points(0, bp[xi],  pch=23)
}

#make a histogram of where the light beams will cross the y axis 
bps = bp
bps[bps<min(fx)] = min(fx)
bp

}


x = seq(-2, 2, by=.01)
xline = seq(-2, 2, by=0.1)
xline[xline==0]=.00001

a = .5
fx = a*x^2  #a parabola
fxline = a*xline^2
fpxline = 2*a*xline

#reflectSurface(x, xline, fx, fxline, fpxline)
a = 1
fx = a*cosh(x/a) #a catenary
fxline = a*cosh(xline/a)
fpxline = sinh(xline/a)

# reflectSurface(x, xline, fx, fxline, fpxline)

a=-1
fx = a*cos(a*x) #a cosine
fxline = a*cos(a*xline)
fpxline = -a^2*sin(a*xline)
#reflectSurface(x, xline, fx, fxline, fpxline)

a=-1
fx = a*sqrt(1-x^2)
fxline = a*sqrt(1-xline^2)
fpxline = -1*a*xline/sqrt(1-xline^2)
reflectSurface(x, xline, fx, fxline, fpxline)