Empirical observations have shown that the product of normal state conductivity and critical temperature is proportional to the superfluid density, which in turn is inversely proportional to the squared penetration depth $\varrho_{s} \propto \frac{1 {\lambda^{2}}$\[1\]. This relation -- known as Homes' scaling -- seems to hold for many conventional and unconventional superconductors. However, its validity and theoretical explanation are still under debate. Furthermore, the superfluid density correlates to the depairing current density $J_{d}$ $\propto \frac{\rho_{s}}{\xi}$ and the core energy of a flux vortex and consequently the available energy gain for efficient flux pinning enabling high critical current densities. Measurements of the critical temperature and normal state conductivity are quite straightforward; however, measurements of the penetration depth are not as trivial, but are highly desirable and may help to elucidate fundamental mechanisms of superconductivity. Different approaches have been reported in literature to measure the penetration depth, ranging from muon spin resonance, mutual inductance measurements, microwave resonance spectroscopy, to magnetization measurements i.e. measuring the magnetic moment of the sample and comparing it to a perfect diamagnetic response with χ=-1. The mathematics is straightforward for an infinite geometry; however, complications arise for finite sample geometries due to demagnetization effects. We present a novel approach for numerically calculating the response of superconducting square thin films in the Meissner phase in a homogenous magnetic field utilizing the fluxoid quantization and taking advantage of symmetries while spanning the grid over our sample. Advantages, limitations, numeric details, the influence of the film thickness and overall geometry will be discussed and put into context of ac SQUID magnetometry experiments.