I can't read those PDFs because you linked to files on your local filesystem, not on the internet.

One way to figure out the problem is with the divide and conquer method. Make sure each part is doing what you expect it to do before you combine them. Maybe even simplify a little.

Some examples: make sure that the Positional_Matrix has the values you expect. What you have here is the same value in every column with no variance in rows, so you might be able to simplify it to be an array until you're sure that the math is right (whether that would help in this context is questionable, but a 2D graph is often easier to understand then a 3D one).
Make sure that the angle and distance orbital calculations are varying smoothly. It's a little harder with an elliptical orbit than a circular one, but not too much.
Make sure that the angle of the sun at each latitude for each season/month varies smoothly from pole to pole (roughly cos(axial_tilt-latitude) to -cos(axial_tilt+latitude) where axial_tilt is the tilt angle pointed toward the sun).

Only after you've gotten to the point of having your insolation at the top of the atmosphere calculation working should you move on to trying to get things like surface heating.