- Timestamp:
- Aug 26, 2010, 9:18:39 AM (16 years ago)
- Location:
- branches/sc_branches/trunkTest
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
ppSim/src/ppSimInsertGalaxies.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/sc_branches/trunkTest
- Property svn:mergeinfo changed
-
branches/sc_branches/trunkTest/ppSim/src/ppSimInsertGalaxies.c
r27657 r29060 1 1 # include "ppSim.h" 2 2 static char *defaultModel = "PS_MODEL_SERSIC"; 3 4 float imageSum (psImage *image) { 5 float sum = 0.0; 6 for (int iy = 0; iy < image->numRows; iy++) { 7 for (int ix = 0; ix < image->numCols; ix++) { 8 sum += image->data.F32[iy][ix]; 9 } 10 } 11 return sum; 12 } 3 13 4 14 bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config) { … … 84 94 psArray *sources = sources = detections->allSources; 85 95 96 // output filename 97 char outname[1024]; 98 char *outroot = psMetadataLookupStr(&mdok, config->arguments, "OUTPUT"); 99 sprintf (outname, "%s.dat", outroot); 100 101 FILE *outfile = fopen (outname, "a"); 102 86 103 // add sources to the readout image & weight 87 104 for (long i = 0; i < galaxies->n; i++) { … … 131 148 132 149 // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) 133 float radius = model->modelRadius (model->params, roughNoise);150 float radius = model->modelRadius (model->params, 0.1*roughNoise); 134 151 radius = PS_MAX (radius, 1.0); 135 152 // XXX the exp(-r^0.25) models can go way out if allowed... 136 radius = PS_MIN (radius, 150.0);153 radius = PS_MIN (radius, 300.0); 137 154 138 155 // construct a source, with model flux pixels set, based on the model 139 156 pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_EXTENDED); 157 158 galaxy->flux = model->modelFlux (model->params); 140 159 141 160 // XXX set the mag & err values (should this be done in pmSourceFromModel?) … … 145 164 source->errMag = sqrt(Area*PS_SQR(roughNoise) + galaxy->flux) / galaxy->flux; 146 165 147 // XXX add the sources to a source array 166 // insert the source flux in the image 167 float sum1 = imageSum(source->pixels); 168 pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); 169 float sum2 = imageSum(source->pixels); 170 float flux = sum2 - sum1; 148 171 149 # if (0) 150 if (CONVOLVED_FIT) { 172 // fprintf (stderr, "flux: %f, sum1: %f, sum2: %f, flux2: %f, dflux: %f\n", galaxy->flux, sum1, sum2, sum2 - sum1, galaxy->flux - sum2 + sum1); 151 173 152 // select the PSF (XXX : move out of loop) 153 pmPSF *psfModel = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF"); 154 assert (psf); 155 156 // supply the psf to the source (normalized?) 157 source->modelPSF = pmModelFromPSFforXY (psfModel, xChip, yChip, 1.0); 174 float par8 = (model->params->n == 8) ? model->params->data.F32[7] : 0.0; 175 fprintf (outfile, "%8.3f %8.3f %10.2f %2d %7.3f %5.3f %5.3f %5.3f %5.3f %5.3f\n", galaxy->x, galaxy->y, flux, 1, source->psfMag, source->errMag, axes.major, axes.minor, axes.theta, par8); 158 176 159 // instantiate the psf flux160 pmSourceCachePSF (source, maskVal);161 162 // convert the cached cached psf model for this source to a psKernel163 // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius)164 // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius)165 psKernel *psf = psphotKernelFromPSF (source, psfSize);166 167 // instantiate the source model flux168 // XXX this may need to handle the offset?169 pmSourceCacheModel (source, maskVal);170 171 // make a storage buffer for the output image172 psImage *buffer = psImageCopy (source->modelFlux, PS_TYPE_F32);173 174 // convolve the psf image with the model image175 psImageConvolveDirect (buffer, source->modelFlux, psf);176 177 // save the result back in the source178 psFree (source->modelFlux);179 source->modelFlux = buffer;180 }181 # endif182 183 // insert the source flux in the image184 pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);185 177 psArrayAdd (sources, 100,source); 186 178 } 179 fclose (outfile); 187 180 188 181 // XXX many leaks in here, i think … … 191 184 return true; 192 185 } 193
Note:
See TracChangeset
for help on using the changeset viewer.
