Polarization of radio pulsar emission contains valuable constraints on the viewing geometry, crucial to interpreting the high energy emission. Classic rotating vector model (RVM) fits often give useful values of magnetic impact angle, beta, for young pulsars. However, RVM fitting generally fails for millisecond pulsars (MSP), where the magnetospheres are small and finite altitude effects are crucial. We have developed a numerical model incorporating multi-component, multi-altitude emission, along with interstellar scattering and emission mode jumps. We show viable fits for several Fermi-detected MSP; the resulting geometry constraints help understand the gamma-ray light curve.