how can i write a program that computes the far field.....

The basic idea is simple enough: convolve the currents with far-field approximation of the Green's function. For multilayer microstrip, you will need the far field approximation of the layered media Green's function. This can be found in Mosig's chapter "Integral Equation Techniques" in "Numerical Techniques for Microwave and Millimetter-Wave Passive Structures", Itoh ed., Wiley-Interscience. If you assume there are no surface waves, it's not to bad. Otherwise, I hope you're familiar with steepest-descent and saddle-point methods (see Felson & Marcwitz, "Radiation and Scattering of Waves," IEEE Press).

I would like to ask taf23 a question that, when you deal with RWG MoM, how do you perform the double integral to get the Z matrix, especially in the case singularities. Do you use Garlerkin method, to "exactly" get the solution, or use the approximate method as provided in the RWG article ?