Changeset 17396 for trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c
- Timestamp:
- Apr 8, 2008, 8:36:06 AM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c
r16819 r17396 3 3 * @author EAM, IfA 4 4 * 5 * @version $Revision: 1. 9$ $Name: not supported by cvs2svn $6 * @date $Date: 2008-0 3-05 01:08:08 $5 * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2008-04-08 18:35:38 $ 7 7 * 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 251 251 } 252 252 253 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname )253 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) 254 254 { 255 255 256 bool status; 256 257 psArray *table; 257 258 psMetadata *row; 258 int i;259 259 psF32 *PAR, *dPAR; 260 psEllipseAxes axes;261 260 psF32 xPos, yPos; 262 261 psF32 xErr, yErr; … … 273 272 table = psArrayAllocEmpty (sources->n); 274 273 274 // which extended source analyses should we perform? 275 bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); 276 bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); 277 bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); 278 bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); 279 280 psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); 281 psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); 282 assert (radialBinsLower->n == radialBinsUpper->n); 283 275 284 // we write out all sources, regardless of quality. the source flags tell us the state 276 for (i = 0; i < sources->n; i++) {285 for (int i = 0; i < sources->n; i++) { 277 286 // skip source if it is not a ext sourc 278 287 // XXX we have two places that extended source parameters are measured: … … 283 292 pmSource *source = sources->data[i]; 284 293 285 // choose the convolved EXT model, if available, otherwise the simple one 286 pmModel *model = source->modelConv; 287 if (model == NULL) { 288 model = source->modelEXT; 289 } 290 if (model == NULL) continue; 291 292 // XXX do we need to do this? 294 // skip sources without measurements 293 295 if (source->extpars == NULL) continue; 296 297 // we require a PSF model fit (ignore the real crud) 298 pmModel *model = source->modelPSF; 299 if (model == NULL) continue; 294 300 295 301 // XXX I need to split the extended models from the extended aperture measurements … … 301 307 yErr = dPAR[PM_PAR_YPOS]; 302 308 303 // XXX for the aperture values, I probably can / should just refer to the psf position and not write any of this304 axes = pmPSF_ModelToAxes (PAR, 20.0);305 306 309 row = psMetadataAlloc (); 310 307 311 // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) 308 312 psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); … … 311 315 psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); 312 316 psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); 313 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag); 314 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); 315 316 // XXX these should be major and minor, not 'x' and 'y' 317 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); 318 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); 319 psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); 320 321 // XXX at the moment, this will be a fixed sized table, but perhaps have multiple output formats depending on what is measured? 322 323 // other values that I need to report: 324 // R, Mag Petrosian, ..... 325 if (source->extpars) { 326 327 // Petrosian measurements 317 318 // Petrosian measurements 319 // XXX insert header data: petrosian ref radius, flux ratio 320 if (doPetrosian) { 328 321 pmSourcePetrosianValues *petrosian = source->extpars->petrosian; 329 322 if (petrosian) { … … 332 325 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); 333 326 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); 327 } else { 328 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); 329 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); 330 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); 331 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); 334 332 } 335 336 // Kron measurements 333 } 334 335 // Kron measurements 336 if (doKron) { 337 337 pmSourceKronValues *kron = source->extpars->kron; 338 338 if (kron) { … … 341 341 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); 342 342 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); 343 } else { 344 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); 345 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); 346 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); 347 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); 343 348 } 344 345 // Isophot measurements 346 // XXX insert header data: isophotal level 349 } 350 351 // Isophot measurements 352 // XXX insert header data: isophotal level 353 if (doIsophotal) { 347 354 pmSourceIsophotalValues *isophot = source->extpars->isophot; 348 355 if (isophot) { … … 351 358 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); 352 359 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); 360 } else { 361 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); 362 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); 363 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); 364 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); 353 365 } 354 366 } 367 368 // Flux Annuli 369 if (doAnnuli) { 355 370 pmSourceAnnuli *annuli = source->extpars->annuli; 356 371 if (annuli) { … … 367 382 sprintf (name, "FLUX_VAR_R_%02d", j); 368 383 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); 369 } 384 } 385 } else { 386 for (int j = 0; j < radialBinsLower->n; j++) { 387 char name[32]; 388 sprintf (name, "FLUX_VAL_R_%02d", j); 389 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); 390 sprintf (name, "FLUX_ERR_R_%02d", j); 391 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); 392 sprintf (name, "FLUX_VAR_R_%02d", j); 393 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); 394 } 370 395 } 371 396 } 372 397 373 374 398 psArrayAdd (table, 100, row); 399 psFree (row); 375 400 } 376 401 377 402 if (table->n == 0) { 378 psFitsWriteBlank (fits, outhead, extname); 379 psFree (table); 380 return true; 403 psFitsWriteBlank (fits, outhead, extname); 404 psFree (outhead); 405 psFree (table); 406 return true; 381 407 } 382 408 383 409 psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); 384 410 if (!psFitsWriteTable (fits, outhead, table, extname)) { 385 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 386 psFree(table); 387 return false; 388 } 411 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 412 psFree (outhead); 413 psFree(table); 414 return false; 415 } 416 psFree (outhead); 389 417 psFree (table); 390 418 … … 397 425 psArray *table; 398 426 psMetadata *row; 399 int i;400 427 psF32 *PAR, *dPAR; 401 428 psEllipseAxes axes; 402 429 psF32 xPos, yPos; 403 430 psF32 xErr, yErr; 431 char name[64]; 404 432 405 433 // create a header to hold the output data … … 412 440 sources = psArraySort (sources, pmSourceSortBySN); 413 441 442 // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams 443 int nParamMax = 0; 444 for (int i = 0; i < sources->n; i++) { 445 pmSource *source = sources->data[i]; 446 if (source->modelFits == NULL) continue; 447 for (int j = 0; j < source->modelFits->n; j++) { 448 pmModel *model = source->modelFits->data[j]; 449 assert (model); 450 nParamMax = PS_MAX (nParamMax, model->params->n); 451 } 452 } 453 414 454 table = psArrayAllocEmpty (sources->n); 415 455 416 456 // we write out all sources, regardless of quality. the source flags tell us the state 417 for (i = 0; i < sources->n; i++) { 418 // skip source if it is not a ext sourc 419 // XXX we have two places that extended source parameters are measured: 420 // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, 421 // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) 422 // should we require both? 457 for (int i = 0; i < sources->n; i++) { 423 458 424 459 pmSource *source = sources->data[i]; 425 460 426 // XXX need to have an array of model fits; loop over the array and write out all 427 // which we calculated 428 429 // choose the convolved EXT model, if available, otherwise the simple one 430 // XXX should not need to choose: write both out 431 pmModel *model = source->modelConv; 432 if (model == NULL) { 433 model = source->modelEXT; 461 // XXX if no model fits are saved, write out modelEXT? 462 if (source->modelFits == NULL) continue; 463 464 // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) 465 for (int j = 0; j < source->modelFits->n; j++) { 466 467 // choose the convolved EXT model, if available, otherwise the simple one 468 pmModel *model = source->modelFits->data[j]; 469 assert (model); 470 471 PAR = model->params->data.F32; 472 dPAR = model->dparams->data.F32; 473 xPos = PAR[PM_PAR_XPOS]; 474 yPos = PAR[PM_PAR_YPOS]; 475 xErr = dPAR[PM_PAR_XPOS]; 476 yErr = dPAR[PM_PAR_YPOS]; 477 478 axes = pmPSF_ModelToAxes (PAR, 20.0); 479 480 row = psMetadataAlloc (); 481 482 // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) 483 psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); 484 psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); 485 psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); 486 psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); 487 psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); 488 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); 489 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); 490 491 psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); 492 psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); 493 494 // XXX these should be major and minor, not 'x' and 'y' 495 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); 496 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); 497 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); 498 499 // write out the other generic parameters 500 for (int k = 0; k < nParamMax; k++) { 501 if (k == PM_PAR_I0) continue; 502 if (k == PM_PAR_SKY) continue; 503 if (k == PM_PAR_XPOS) continue; 504 if (k == PM_PAR_YPOS) continue; 505 if (k == PM_PAR_SXX) continue; 506 if (k == PM_PAR_SXY) continue; 507 if (k == PM_PAR_SYY) continue; 508 509 snprintf (name, 64, "EXT_PAR_%02d", k); 510 511 if (k < model->params->n) { 512 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); 513 } else { 514 psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); 515 } 516 } 517 518 // XXX other parameters which may be set. 519 // XXX flag / value to define the model 520 // XXX write out the model type, fit status flags 521 522 psArrayAdd (table, 100, row); 523 psFree (row); 434 524 } 435 if (model == NULL) continue;436 437 PAR = model->params->data.F32;438 dPAR = model->dparams->data.F32;439 xPos = PAR[PM_PAR_XPOS];440 yPos = PAR[PM_PAR_YPOS];441 xErr = dPAR[PM_PAR_XPOS];442 yErr = dPAR[PM_PAR_YPOS];443 444 axes = pmPSF_ModelToAxes (PAR, 20.0);445 446 row = psMetadataAlloc ();447 // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)448 psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq);449 psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos);450 psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos);451 psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr);452 psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr);453 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag);454 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag);455 456 // XXX these should be major and minor, not 'x' and 'y'457 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major);458 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor);459 psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta);460 461 // XXX other parameters which may be set.462 // XXX flag / value to define the model463 464 psArrayAdd (table, 100, row);465 psFree (row);466 525 } 467 526 468 527 if (table->n == 0) { 469 psFitsWriteBlank (fits, outhead, extname); 470 psFree (table); 471 return true; 528 psFitsWriteBlank (fits, outhead, extname); 529 psFree (outhead); 530 psFree (table); 531 return true; 472 532 } 473 533 474 534 psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); 475 535 if (!psFitsWriteTable (fits, outhead, table, extname)) { 476 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 477 psFree(table); 478 return false; 479 } 536 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 537 psFree (outhead); 538 psFree(table); 539 return false; 540 } 541 psFree (outhead); 480 542 psFree (table); 481 482 543 return true; 483 544 }
Note:
See TracChangeset
for help on using the changeset viewer.
