If the ranges are not too long (so that we can neglect the curvature of the earth) and if altitude differences can be ignored, I would start by converting the lat/lon coords of the sensor at the times of measurement to a local Cartesian x-y system. Do this by taking a lat/lon measurement near the centroid of your sensor measurements. That reference measurement will be the origin of your local x-y system, i.e. (0,0). Then for every other lat/lon measurement compute the range (rng) and bearing (brg) from the ref lat/lon. Then convert the rng/brg coords to to x-y coords.

OK so now we have a 2D problem, with a set of brg measurements at x-y points. For a reasonable estimate at a low computational cost, just use the pseudo-linear estimator (PLE). This yields an estimate of the source location, using a simple least-squares type approach, requiring the inverse of a matrix. For a better estimate, if you have a reasonable idea of the expected error in your bearing measurements, and you don't mind doing a bit of coding, use the maximum likelihood estimator (MLE). This is an iterative approach solved using the Gauss-Newton procedure. Initialize it using the PLE estimate.

Spingarn K. Passive position location estimation using the extended Kalman filter. IEEE Trans Aerosp Electron Syst AES 1987;23:558–67.

It is not too difficult to put together some C code to do all of this, but I get the feeling you are looking for a ready-to-go tool.