When near the focus, the image of a point emitter has little information about the emitter's position along the optical axis. To improve the 3D localization precision, several methods such as astigmatic imaging, double helix point spread functions, and multiple focal plane strategies have been developed to encode the axial position in one or more 2D images of the same emitter. In each of these strategies, a proper model for the 3D microscope point spread function (PSF) is essential to avoid inaccuracies in localization. Generating PSF models from phase retrieved (PR) pupil functions has the advantage that the model images are noise free, can be calculated for any focal plane, and include aberrations specific to the optical system. This representation can also be extremely compact by representing the pupil function as a set of Zernike polynomials. Here we demonstrate the use of PSFs generated from PR pupil functions for 3D single molecule localization based super-resolution imaging using a dual focal plane setup. We show the phase retrieval process, as well as the details of fast PSF calculation by implementation on graphics processing unit hardware. We show a comparison between fitting with ideal, aberration free PSF models and that generated from PR pupil functions using a variety of cellular structures including the ER, microtubules and the membrane proteins STIM1 and ORAI1.