This program works well for the example given but when I plot a NACA 65016 their is an obvious discontinuity at the point of max camber. I believe the equations used are correct, but I'm not exactly sure how the k1 constant is obtained. Either way I believe the z-values at the point of max camber are wrong.
