#include #include #include #include #include #include "fitsio.h" #include "longnam.h" #include "structs.h" void objects_head (struct fitpars *fpar, fitsfile *fptr, int offx, int offy); void printerror(int status); int istrunc (char *objtype); void writeheader (struct image *img, struct fitpars *fpar, double chisq, int ndof, int nfree, int nfixed, int hdu, int offx, int offy) { extern int (*pfunc) (const char *, ...); extern struct inpars input; char param[FLEN_KEYWORD]; fitsfile *fptr; int status = 0; float chi2nu; if (fits_open_file (&fptr, img->name, READWRITE, &status)) { pfunc ("\n Can't append to image header %s!\n", img->name); printerror( status ); }; if (fits_movabs_hdu (fptr, hdu, IMAGE_HDU, &status)) { pfunc ("\n Can't append to image header %s!\n", img->name); printerror( status ); }; fits_write_comment (fptr, "", &status); fits_write_comment (fptr, "========== GALFIT Input Parameters ==========", &status); fits_write_comment (fptr, "", &status); fits_update_key(fptr, TSTRING, "INITFILE", input.initparfile, "GALFIT input file", &status); fits_update_key(fptr, TSTRING, "DATAIN", input.data, "Input data image", &status); fits_update_key(fptr, TSTRING, "SIGMA", input.sigma, "Input sigma image", &status); sprintf (param, "%s %s", input.psf, input.kernel); fits_update_key(fptr, TSTRING, "PSF", param, "Convolution PSF and kernel", &status); fits_update_key(fptr, TSTRING, "CONSTRNT", input.constraints, "Parameter constraint file", &status); fits_update_key(fptr, TSTRING, "MASK", input.badpix, "Input mask image", &status); fits_update_key(fptr, TSTRING, "FITSECT", input.imgsect, "Image section fitted", &status); sprintf (param, "%i, %i", input.convbox[1], input.convbox[2]); fits_update_key(fptr, TSTRING, "CONVBOX", param, "Convolution box size", &status); fits_update_key(fptr, TFLOAT, "MAGZPT", &input.magzpt, "Magnitude zeropoint", &status); fits_write_comment (fptr, "========== GALFIT Final Parameters ==========", &status); objects_head (fpar, fptr, offx, offy); if (!isnormal((float)chisq)) chisq = 0.; fits_update_key(fptr, TDOUBLE, "Chisq", &chisq, "Chi^2 of fit", &status); fits_update_key(fptr, TINT, "NDOF", &ndof, "Degrees of freedom", &status); fits_update_key(fptr, TINT, "NFREE", &nfree, "Number of free parameters", &status); fits_update_key(fptr, TINT, "NFIX", &nfixed, "Number of fixed parameters", &status); if (ndof > 0) chi2nu = chisq/ndof; else chi2nu = 0.; fits_update_key(fptr, TFLOAT, "Chi2nu", &chi2nu, "Reduced Chi^2", &status); fits_update_key(fptr, TSTRING, "LOGFILE", img->logfile, "Output logfile", &status); fits_write_comment (fptr, "=============================================", &status); fits_write_comment (fptr, " ", &status); fits_close_file(fptr, &status); /* close the file */ return; } void objects_head (struct fitpars *fpar, fitsfile *fptr, int offx, int offy) { char par[FLEN_KEYWORD], skyc[FLEN_KEYWORD], comp[FLEN_KEYWORD], key[FLEN_KEYWORD], comment[FLEN_COMMENT]; float parval; struct fitpars *hptr; extern float xskycent, yskycent; int status = 0, i, j, k; j = 1; while (fpar != NULL) { sprintf (comp, "COMP_%i", j); fits_update_key(fptr, TSTRING, comp, fpar->objtype, "Component type", &status); /***********************************\ * First do classical parameters * \***********************************/ for (i=0; i<= 10; i++) { if (i==1 && strncmp(fpar->objtype, "sky", 3) != 0) { if (fpar->ia[1] != -1) parval = fpar->a[1] + offx; } else if (i==2 && strncmp(fpar->objtype, "sky", 3) != 0) { if (fpar->ia[2] != -1) parval = fpar->a[2] + offy; } else if (i != 0) parval = fpar->a[i]; if (i != 0 && fabs(parval) > 0.01) { if (fpar->ia[i] == 1) sprintf (par, "%-.4f +/- %-.4f", parval, fabs(fpar->sig[i])); else if (fpar->ia[i] == 0) sprintf (par, "[%-.4f]", parval); else if (fpar->ia[i] == 2) sprintf (par, "{%-.4f}", parval); } else if (i != 0) { if (fpar->ia[i] == 1) sprintf (par, "%-.3e +/- %-.3e", parval, fabs(fpar->sig[i])); else if (fpar->ia[i] == 0) sprintf (par, "[%-.3e]", parval); else if (fpar->ia[i] == 2) sprintf (par, "{%-.3e}", parval); }; if (strncmp(fpar->objtype, "sky", 3) != 0) { if (i==1 && (fpar->ia[1] >= 0 || fpar->ia[2] >= 0)) { sprintf (key, "%i_XC", j); fits_update_key(fptr, TSTRING, key, par, "X center [pixel]", &status); }; if (i==2 && (fpar->ia[1] >= 0 || fpar->ia[2] >= 0)) { sprintf (key, "%i_YC", j); fits_update_key(fptr, TSTRING, key, par, "Y center [pixel]", &status); }; /*******\ * Nuker * \*******/ if (strncmp(fpar->objtype, "nuker", 5) == 0) { if (i==3) { if (index(fpar->objtype, '/') != NULL) { sprintf (key, "%i_MUBRK", j); fits_update_key(fptr, TSTRING, key, par, "Surf. brghtnss @ break radius [mag/arcsec^2]", &status); } else if (index(fpar->objtype, '\\') != NULL ) { sprintf (key, "%i_MU_0", j); fits_update_key(fptr, TSTRING, key, par, "Central surface brightness [mag/arcsec^2]", &status); } else { sprintf (key, "%i_MU", j); fits_update_key(fptr, TSTRING, key, par, "Surface brightness at Rb", &status); }; }; if (i==4) { sprintf (key, "%i_Rb", j); fits_update_key(fptr, TSTRING, key, par, "Break radius Rb [pixels]", &status); }; if (i==5) { sprintf (key, "%i_alpha", j); fits_update_key(fptr, TSTRING, key, par, "alpha", &status); }; if (i==6) { sprintf (key, "%i_beta", j); fits_update_key(fptr, TSTRING, key, par, "beta", &status); }; if (i==7) { sprintf (key, "%i_gamma", j); fits_update_key(fptr, TSTRING, key, par, "gamma", &status); }; }; /**************************\ * Sersic and deVaucouleurs * \**************************/ if (strncmp(fpar->objtype, "sersic", 6) == 0 || strncmp(fpar->objtype, "devauc", 6) == 0) { if (i==3) { if (index(fpar->objtype, '/') != NULL) { sprintf (key, "%i_MUBRK", j); fits_update_key(fptr, TSTRING, key, par, "Surf. brghtnss @ break radius [mag/arcsec^2]", &status); } else if (index(fpar->objtype, '\\') != NULL ) { sprintf (key, "%i_MU_0", j); fits_update_key(fptr, TSTRING, key, par, "Central surface brightness [mag/arcsec^2]", &status); } else if (strncmp(fpar->objtype, "sersic2", 7) == 0 || strncmp(fpar->objtype, "devauc2", 7) == 0) { sprintf (key, "%i_MU_E", j); fits_update_key(fptr, TSTRING, key, par, "Effective surface brightness [mag/arcsec^2]", &status); } else { sprintf (key, "%i_MAG", j); fits_update_key(fptr, TSTRING, key, par, "Integrated magnitude [mag]", &status); }; }; if (i==4) { sprintf (key, "%i_Re", j); fits_update_key(fptr, TSTRING, key, par, "Effective radius Re [pixels]", &status); }; if (i==5 && strncmp(fpar->objtype, "sersic", 6) == 0) { sprintf (key, "%i_n", j); fits_update_key(fptr, TSTRING, key, par, "Sersic index", &status); }; }; /******************\ * Exponential Disk * \******************/ if (strncmp(fpar->objtype, "expdisk", 7) == 0) { if (i==3) { if (index(fpar->objtype, '/') != NULL) { sprintf (key, "%i_MUBRK", j); fits_update_key(fptr, TSTRING, key, par, "Surf. brghtnss @ the break radius [mag/arcsec^2]", &status); } else if (index(fpar->objtype, '\\') != NULL ) { sprintf (key, "%i_MU_0", j); fits_update_key(fptr, TSTRING, key, par, "Central surface brightness [mag/arcsec^2]", &status); } else { sprintf (key, "%i_MAG", j); fits_update_key(fptr, TSTRING, key, par, "Integrated magnitude", &status); }; }; if (i==4) { sprintf (key, "%i_Rs", j); fits_update_key(fptr, TSTRING, key, par, "Scalelength [pixels]", &status); }; }; /**************\ * Edge-on Disk * \**************/ if (strncmp(fpar->objtype, "edgedisk", 8) == 0) { if (i==3) { sprintf (key, "%i_MU_0", j); fits_update_key(fptr, TSTRING, key, par, "Central surface brightness", &status); }; if (i==4) { sprintf (key, "%i_Hs", j); fits_update_key(fptr, TSTRING, key, par, "Scale-height [pixels]", &status); }; if (i==5) { sprintf (key, "%i_Rs", j); fits_update_key(fptr, TSTRING, key, par, "Scale-length [pixels]", &status); }; }; /********\ * Ferrer * \********/ if (strncmp(fpar->objtype, "ferrer", 6) == 0) { if (i==3) { sprintf (key, "%i_MU_0", j); fits_update_key(fptr, TSTRING, key, par, "Central surface brightness", &status); }; if (i==4) { sprintf (key, "%i_Rout", j); fits_update_key(fptr, TSTRING, key, par, "Outer radius [pixels]", &status); }; if (i==5) { sprintf (key, "%i_ALPHA", j); fits_update_key(fptr, TSTRING, key, par, "Alpha (outer truncation sharpness)", &status); }; if (i==6) { sprintf (key, "%i_BETA", j); fits_update_key(fptr, TSTRING, key, par, "Beta (central slope)", &status); }; }; /*********************\ * Gaussian and Moffat * \*********************/ if (strncmp(fpar->objtype, "gauss", 5) == 0 || strncmp(fpar->objtype, "moffat", 6) == 0 || strncmp(fpar->objtype, "psf", 3) == 0 ) { if (i==3) { if (index(fpar->objtype, '/') != NULL) { sprintf (key, "%i_MUBRK", j); fits_update_key(fptr, TSTRING, key, par, "Surf. brghtnss @ the break radius [mag/arcsec^2]", &status); } else if (index(fpar->objtype, '\\') != NULL ) { sprintf (key, "%i_MU_0", j); fits_update_key(fptr, TSTRING, key, par, "Central surface brightness [mag/arcsec^2]", &status); } else { sprintf (key, "%i_MAG", j); fits_update_key(fptr, TSTRING, key, par, "Integrated magnitude", &status); }; }; if (i==4 && strncmp(fpar->objtype, "psf", 3) != 0) { sprintf (key, "%i_FWHM", j); fits_update_key(fptr, TSTRING, key, par, "FWHM [pixels]", &status); }; if (i==5 && strncmp(fpar->objtype, "moffat", 6) == 0) { sprintf (key, "%i_C", j); fits_update_key(fptr, TSTRING, key, par, "Powerlaw", &status); }; }; /******\ * King * \******/ if (strncmp(fpar->objtype, "king", 4) == 0) { if (i==3) { if (index(fpar->objtype, '/') != NULL) { sprintf (key, "%i_MUBRK", j); fits_update_key(fptr, TSTRING, key, par, "Surf. brghtnss @ the break radius [mag/arcsec^2]", &status); } else if (index(fpar->objtype, '\\') != NULL ) { sprintf (key, "%i_MU_0", j); fits_update_key(fptr, TSTRING, key, par, "Central surface brightness [mag/arcsec^2]", &status); }; }; if (i==4) { sprintf (key, "%i_Rc", j); fits_update_key(fptr, TSTRING, key, par, "Core radius [pixels]", &status); }; if (i==5) { sprintf (key, "%i_Rt", j); fits_update_key(fptr, TSTRING, key, par, "Truncation radius [pixels]", &status); }; if (i==6) { sprintf (key, "%i_alpha", j); fits_update_key(fptr, TSTRING, key, par, "Powerlaw alpha", &status); }; }; }; /************\ * Truncation * \************/ if (istrunc(fpar->objtype)) { if (i == 0 && fpar->inner != NULL) { sprintf (key, "%i_INNER", j); sprintf (par, ""); for (k = 1; k <= fpar->inner[0]; k++) sprintf (par, "%s %d", par, fpar->inner[k]); if (fpar->inner[0] == 0) sprintf (par, "0"); fits_update_key(fptr, TSTRING, key, par, "Component(s) with inner truncation", &status); }; if (i == 0 && fpar->outer != NULL) { sprintf (key, "%i_OUTER", j); sprintf (par, ""); for (k = 1; k <= fpar->outer[0]; k++) sprintf (par, "%s %d", par, fpar->outer[k]); if (fpar->outer[0] == 0) sprintf (par, "0"); fits_update_key(fptr, TSTRING, key, par, "Component(s) with outer truncation", &status); }; if (i==4) { sprintf (key, "%i_RBRK", j); fits_update_key(fptr, TSTRING, key, par, "Break radius (99\% normal flux) [pixels]", &status); }; if (i==5) { if (index (fpar->objtype, '2') != NULL) { sprintf (key, "%i_RSFT", j); fits_update_key(fptr, TSTRING, key, par, "Softening radius (1\% normal flux) [pixels]", &status); } else { sprintf (key, "%i_DSFT", j); fits_update_key(fptr, TSTRING, key, par, "Softening length (1\% normal flux) [pixels]", &status); }; }; }; /*******************\ * Common parameters * \*******************/ if (i==9 && strncmp(fpar->objtype, "psf", 3) != 0 && strncmp(fpar->objtype, "edgedisk", 8) != 0 && fpar->ia[9] != -1) { sprintf (key, "%i_AR", j); fits_update_key(fptr, TSTRING, key, par, "Axis ratio (b/a)", &status); }; if (i==10 && strncmp(fpar->objtype, "psf", 3) != 0 && fpar->ia[10] != -1) { sprintf (key, "%i_PA", j); fits_update_key(fptr, TSTRING, key, par, "Position Angle (PA) [degrees: Up=0, Left=90]", &status); }; /*****\ * Sky * \*****/ if (strncmp(fpar->objtype, "sky", 3) == 0) { if (i==1) { sprintf (skyc, "[%-.4f]", xskycent + offx); sprintf (key, "%i_XC", j); fits_update_key(fptr, TSTRING, key, skyc, "X center [pixel]", &status); sprintf (skyc, "[%-.4f]", yskycent + offy); sprintf (key, "%i_YC", j); fits_update_key(fptr, TSTRING, key, skyc, "Y center [pixel]", &status); sprintf (key, "%i_SKY", j); fits_update_key(fptr, TSTRING, key, par, "Sky background [ADUs]", &status); }; if (i==2) { sprintf (key, "%i_DSDX", j); fits_update_key(fptr, TSTRING, key, par, "x sky gradient [ADUs]", &status); }; if (i==3) { sprintf (key, "%i_DSDY", j); fits_update_key(fptr, TSTRING, key, par, "y sky gradient [ADUs]", &status); }; }; }; /**************************************\ * Truncation flag for light profiles * \**************************************/ if (!istrunc(fpar->objtype)) { if (fpar->inner != NULL) { sprintf (par, ""); for (k = 1; k <= fpar->inner[0]; k++) sprintf (par, "%s %d ", par, fpar->inner[k]); sprintf (key, "%i_TINNR", j); fits_update_key(fptr, TSTRING, key, par, "Inner truncation by component(s)", &status); }; if (fpar->outer != NULL) { sprintf (par, ""); for (k = 1; k <= fpar->outer[0]; k++) sprintf (par, "%s %d ", par, fpar->outer[k]); sprintf (key, "%i_TOTR", j); fits_update_key(fptr, TSTRING, key, par, "Outer truncation by component(s)", &status); }; }; /***************************\ * Higher order parameters * \***************************/ if (fpar->high != NULL) { hptr = fpar->high; while (hptr != NULL) { /**********************\ * Diskyness/Boxyness * \**********************/ if (strncmp (hptr->objtype, "C0", 2) == 0) { parval = hptr->a[1]; if (hptr->ia[1] == 1) sprintf (par, "%-.4f +/- %-.4f", parval, hptr->sig[1]); else if (hptr->ia[1] == 0) sprintf (par, "[%-.4f]", parval); else if (hptr->ia[1] == 2) sprintf (par, "{%-.4f}", parval); sprintf (key, "%i_C0", j); fits_update_key(fptr, TSTRING, key, par, "Diskyness (<0) or Boxyness (>0)", &status); }; /*****************\ * Bending modes * \*****************/ if (strncmp (hptr->objtype, "bend", 4) == 0) { for (i=1; i <= hptr->npars; i++) { if (hptr->ia[i] > -1) { parval = hptr->a[i]; sprintf (comment, "Bending mode %d amplitude", i); sprintf (key, "%i_B%d", j, i); if (hptr->ia[i] == 1) sprintf (par, "%-.4f +/- %-.4f", parval, hptr->sig[i]); else if (hptr->ia[i] == 0) sprintf (par, "[%-.4f]", parval); else if (hptr->ia[i] == 2) sprintf (par, "{%-.4f}", parval); fits_update_key(fptr, TSTRING, key, par, comment, &status); }; }; }; /*****************\ * Fourier modes * \*****************/ if (strncmp (hptr->objtype, "four", 4) == 0) { for (i=1; i <= hptr->npars; i+=2) { if (hptr->ia[i] > -1) { parval = hptr->a[i]; sprintf (comment, "Azimuthal Fourier mode %d amplitude", i/2+1); sprintf (key, "%i_F%d", j, i/2+1); if (hptr->ia[i] == 1) sprintf (par, "%-.4f +/- %-.4f", parval, hptr->sig[i]); else if (hptr->ia[i] == 0) sprintf (par, "[%-.4f]", parval); else if (hptr->ia[i] == 2) sprintf (par, "{%-.4f}", parval); fits_update_key(fptr, TSTRING, key, par, comment, &status); parval = hptr->a[i+1]; sprintf (comment, "Azimuthal Fourier mode %d phase angle", i/2+1); sprintf (key, "%i_F%dpa", j, i/2+1); if (hptr->ia[i+1] == 1) sprintf (par, "%-.4f +/- %-.4f", parval, hptr->sig[i+1]); else if (hptr->ia[i+1] == 0) sprintf (par, "[%-.4f]", parval); else if (hptr->ia[i+1] == 2) sprintf (par, "{%-.4f}", parval); fits_update_key(fptr, TSTRING, key, par, comment, &status); }; }; }; /**********************\ * Rotation functions * \**********************/ if (strncmp (hptr->objtype, "power", 5) == 0 || strncmp (hptr->objtype, "log", 3) == 0) { sprintf (key, "%i_ROTF", j); sprintf (par, "%s", hptr->objtype); fits_update_key(fptr, TSTRING, key, par, "Axial rotation function", &status); for (i=1; i <= hptr->npars; i++) { if (hptr->ia[i] == 1) sprintf (par, "%-.4f +/- %-.4f", hptr->a[i], hptr->sig[i]); else if (hptr->ia[i] == 0) sprintf (par, "[%-.4f]", hptr->a[i]); else if (hptr->ia[1] == 2) sprintf (par, "{%-.4f}", hptr->a[i]); if (i==1) { sprintf (key, "%i_RIN", j); fits_update_key(fptr, TSTRING, key, par, "Inner (const. PA) radius [pixels]", &status); }; if (i==2) { sprintf (key, "%i_ROUT", j); fits_update_key(fptr, TSTRING, key, par, "Outer (const. PA) radius [pixels]", &status); }; if (i==3) { sprintf (key, "%i_RANG", j); fits_update_key(fptr, TSTRING, key, par, "Total rotation angle out to Rout [deg]", &status); }; if (i==4) { if (strncmp (hptr->objtype, "power", 5) == 0) { sprintf (key, "%i_ALPHA", j); fits_update_key(fptr, TSTRING, key, par, "Asymptotic spiral powerlaw", &status); }; if (strncmp (hptr->objtype, "log", 3) == 0) { sprintf (key, "%i_RWS", j); fits_update_key(fptr, TSTRING, key, par, "Logarithmic winding scale radius [pixels]", &status); }; }; if (i==9) { sprintf (key, "%i_INCL", j); fits_update_key(fptr, TSTRING, key, par, "Inclination angle to line of sight [deg]", &status); }; if (i==10) { sprintf (key, "%i_SPA", j); fits_update_key(fptr, TSTRING, key, par, "Projected disk position angle [deg]", &status); }; }; }; hptr = hptr->next; }; }; fits_write_comment (fptr, "------------------------", &status); fpar = fpar->next; j++; }; }