Please do not use this routine. Instead use freduceimage2, which doesn't use any global (common block) variables, and hence is much safer, easier to reuse, and easier to maintain.
Initial setup; if "image" is set, then use freadimage2,/noimage to read the image header information into ims. Otherwise, assume that ims already has that information loaded.
Check the reducedimages database for a NetCDF file of the appropriate information. If it's there, read the list of {star_struct} variables into pho. Return.
If not, and no image data was passed, read and autosurface the image into im. (This will only work for images in the database.)
Run fisofind2 - Loads pho with catalog of objects found.
Run fcatalog2 - Load pho.type with object type and ims.fwx and ims.fwy with a (not too great) estimate for the seeing of the image.
Run faperture2 - Measures the centroid (in pho.x and pho.y), aperture photometry in a ims.fwx x ims.fwy (returned from fcatalog2) aperture (in pho.aper, uncertainty pho.eaper).
If "calcringsig" was set (it is by default), run ringsig2
If tele is 'kp4m', run fcatalogpsp2 - Fits gaussians over things, loads ims.sig, sets several fields in pho.
Do more scary ringsig stuff I don't understand yet.
Run faperture again, with radius=2.36*ims.sig[0],/norec,/usexandy.
Run apmmatch2
-- loads the apm variables in ims.If "nosave" was not set ("nodb" implies "nosave"), then...
Please avoid using fisofind (without the "2") since that uses common block global variables that generally make a huge mess out of things.
NOTE: fisofind2 assumes the image is surfaced!!!
Call the external routine isofind in isofind.so to do an isophotal finding on the image.
If isofind reached it's limit (10,000 stars?), load up ims.comment with a warning to that effect.
create a {star_struct} array, to return in "star", and load it up with the results returned from isofind.so
Clip out and keep found objects with: isophot>mincounts, mxx and mxy != 0, central pixel > 0 (im[round(cmx),round(cmy)]>0).
Also optionally return "worka" which is a lonarr the same size as the input image which tells you which pixels belong to which object. For example, all of the pixels which belonged to object 314 have a value of 314 in worka. This also contains objects which did not pass all of the tests.
NOTE: Assumes surfaced image!!!!
Find pho[i].fwx and pho[i].fwy by searching on each axis, finding the pixel that causes im[x,y] to fall to half the value at the central pixel (im[round(pho[i].cmx),round(pho[i].cmy)]), and linearly interpolating (? -- check this) across that pixel boundary. Take the average in each direction. Set fwx and fwy to 0 if this takes off the edge of the image.
Preliminarily set all images as galaxies (pho.type='gal')
Define cosmic rays where pho.npix<7 (from isofind) and pho.max/ims.exp_sigma>20; set those pho.type to 'cr'
Find the median pho.myy/pho.mxx; find all objects with pho.mxy/pho.isophot<0.2, and pho.myy/pho.mxx within 0.2 of the median of that ration; define ims.fwx and ims.fwy as the median of those values for the objects which match these criterea.
Define as 'stlr' all non-cr objects with pho.fwx and pho.fwy>0 and pho.fwx<1.2*ims.fwx and pho.fwy<1.2*fwy, and pho.mxy/pho.moment<0.2
Define as 'strk' (streaks) all objects with sqrt(pho.moment/pho.isophot*12)>300 (???)
Define as 'sat' all objects with their max pixel value (pho.max, as returned from isofind) > 0.85*ims.welldepth.
Define all other objects as galaxies ('gal').
For telescopes other than 'kp4m', call catalogpsp. Set ims.ringsig to ims.sig[0], ims.errringsig to ims.sig[1], ims.ringrad to [0,1,2,3,4,10000} (?????!?!?!?!???), and ims.nringsig to 6.
For kp4m, do a whole lot of scary stuff that I wish wasn't there, and which is a horrible hack that should really be obviated by Michael Wood-Vasey's current work on a general method of doing subtractions with a psf (or, equivalently, convolution kernel) that varies across an image.
Notes: Uses the "worka" array returned from isofind.
Set variable "rad" to the average of ims.fwx and ims.fwy.
Select appropriate objects in the passed pho array to use: aper>0, aper/eaper>25, type!='sat', type!='cr', npix>5, objects no closer than rad+2 to the edge of the image.
If no objects remain, set ims.sig[0] to zero and return.
If fewer than 50 objects remain and /tryfew is set (freduceimage sets it by default), then try objects with aper/eaper>10, this time including those with pho.type='cr'. (Is this intentional??)
If fewer than 50 objects still remain, set those object's pho.type to 'stand', set ims.sig[*] to 0, and give up.
Fit a 2d parabola over the top of each object in order to calculate pho.alpha, which seems (check this) to be some sort of position angle. (???)
Find a mode of alpha, using the mode algorithm with a window of 25, setting maxmin to the values that fall within that window. Set lastu to the elements of pho that fall within that window. (????).
Iterate up to 5 times, looking for ims[a].sig[0] to converge, on the next several steps
Call buildpsp on the lastu objects
sigma[0] = central valueOverall sigma = sigma[0] + sigma[1]*(cx-nx/2) + sigma[2]*(cy-ny/2)
sigma[1] = x coefficient
sigma[2] = y coefficient
Call fitpsp on all objects remaining in in pho (which is a superset of lastu)
Select the 50 objects with the greatest pspsharp, and of those the 25 with the lowest psperrfit. Set lastu to this selection, and iterate.
Select standard stars, setting pho.type='stand', based on those who pass within a clipped mean about the mode of pho.pspsharp (and, seemingly (oddly????) of pho.psperrfit). (?!?!?!?)
NOTES: I'm a little worried about that clipped mean on pho.pspsharp, since not all elements of pho will have been passed through buildpsp...!