Help with solving an unsteady problem with zero Neumann boundary condition in a circular region.

Hi guys,

I have a PDE/ODE system which I have to solve with finite difference method. The PDE is an unsteady equation which evolves through time. The problem is that the region is circular and in some mesh points I need to have the normal vector (it's not (0,1) nor (1,0)) in order to apply the Neumann boundary condition. I mean is there any algorithms or code which calculates the normal vector in any given points on boundary?

I desperately need your help.