The field in the plane of your detector is the sum of the reference field and the scattered field. The recorded hologram is the intensity of this superposition. Typically, the reference field is easy to describe mathematically. The field scattered by the object, however, can be very complicated.
It would help to know what kind of object you'd like to simulate. If the object is comparable in size to the wavelength of light, (Generalized) Lorenz-Mie theory is useful for calculating the scattered wave. Scattering by more complicated structures can be calculated with the discrete dipole approximation (DDA). If the object is large, then you might be better off using the FDTD (finite-difference time domain) algorithm. Naturally, reflection is different from transmission. In all of these cases, the calculation has to into account the polarization of the illumination, and so yields the scattered wave as a vector function. Fortunately, codes for all of these algorithms are freely available on the web in a variety of languages. Unfortunately, the associated programming tends to be pretty detailed and time-consuming.
You should read in the polarization gratings subjectm, where the holograms is gnerated from the superposition of two plane waves with mutually orthogonal polarizations and equal amplitudes. You could write the whole system using Jones vectors and matrix.
There is a PhD thesis about polarization imaging for digital holographic microscopy. I haven't studied the polarization part in depth, but it may give you some hints about how to simulate it.
Simulating a linearly polarized object beam with a linearly polarized reference beam shouldn't be too difficult. Create a 2D complex field with two components (parallel and orthogonal) for the object beam. Add them together with different ratios to simulate a rotation of the polarization direction on the reference beam. Then add to the 2D complex reference distribution and calculate the intensity.