#include #include #include #include #include "structs.h" int istrunc (char *objtype); float pa (int mode, float instr_pa); int par2num (char *parname); int assign_pars (struct inpars *input, struct image *data, struct image *sigma, struct image *psf, struct image *badpix, struct fitpars *fpar, struct convpars *cpar, struct cons *constr) { int i, mode, parnum, m, itmp, err = 0; float tmp; struct fitpars *ptr, *fparptr; strcpy (data->name, input->data); strcpy (sigma->name, input->sigma); strcpy (psf->name, input->psf); strcpy (badpix->name, input->badpix); if (strncmp (input->kernel, "#", 1) == 0) sprintf (input->kernel, ""); /* code fparptr->ia[i] = -1 means this parameter is not being used */ fparptr = fpar; while (fparptr != NULL) { if (strncmp (fparptr->objtype, "nuker", 5)==0) { fparptr->ia[8] = -1; }; if (strncmp (fparptr->objtype, "sersic", 6)==0) for (i=6; i<=8 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "devauc", 6)==0) { fparptr->a[5] = 4.; fparptr->ia[5] = 0; for (i=6; i<=8 ; i++) fparptr->ia[i] = -1; }; if (strncmp (fparptr->objtype, "isophote", 8)==0) { fparptr->ia[4] = 0; fparptr->ia[5] = 0; for (i=7; i<=8 ; i++) fparptr->ia[i] = -1; }; if (strncmp (fparptr->objtype, "expdisk", 7)==0) for (i=5; i<=8 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "edgedisk", 8)==0) for (i=6; i<=9 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "sinc", 4)==0) for (i=5; i<=8 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "gaussian", 8)==0) for (i=5; i<=8 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "psf", 3)==0) for (i=4; i<=10 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "king", 4)==0) { if (fparptr->a[6] == 0.) { fparptr->a[6] = 2.; fparptr->ia[6] = 0; }; for (i=7; i<=8 ; i++) fparptr->ia[i] = -1; }; if (strncmp (fparptr->objtype, "moffat", 6)==0) for (i=6; i<=8 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "ferrer", 6)==0) for (i=7; i<=8 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "sky", 3)==0) for (i=4; i<=10 ; i++) fparptr->ia[i] = -1; if (strncmp (fparptr->objtype, "power", 5)==0 || strncmp (fparptr->objtype, "log", 3)==0) fparptr->ia[0] = 0; if (istrunc (fparptr->objtype)) { if (index (fparptr->objtype, '2') != NULL) { /*------------------------------------------------------------\ | Type 2 trunc: When the softening free parameter is the | | softening *radius* and instead of softening length. | \------------------------------------------------------------*/ if (fparptr->a[4] > fparptr->a[5]) { tmp = fparptr->a[5]; fparptr->a[5] = fparptr->a[4]; fparptr->a[4] = tmp; itmp = fparptr->ia[5]; fparptr->ia[5] = fparptr->ia[4]; fparptr->ia[4] = itmp; }; }; fparptr->ia[3] = -1; for (i=6; i<=8 ; i++) fparptr->ia[i] = -1; if (fparptr->ia[9] == -1) fparptr->ia[10] = -1; if (fparptr->ia[9] == 1 || fparptr->ia[9] == 0) if (fparptr->ia[10] <= -1 || fparptr->ia[10] > 1) fparptr->ia[10] = 0; }; /* Do something about PA and axis ratio */ if (strncmp (fparptr->objtype, "sky", 3) != 0) { if ( fparptr->a[9] > 1.) { /* Size needs to be multiplied by 1/q */ fparptr->a[4] = fparptr->a[4] * fparptr->a[9]; ptr = fparptr->high; while (ptr != NULL) { if (istrunc (ptr->objtype)) { ptr->a[4] = ptr->a[4] * fparptr->a[9]; ptr->a[5] = ptr->a[5] * fparptr->a[9]; }; ptr = ptr->next; }; /* Invert the axis ratio */ fparptr->a[9] = 1./ fparptr->a[9]; /* Flip the PA by 90 deg. */ fparptr->a[10] = fparptr->a[10] + 90.; }; mode = 2; /* coord rotation has to retain PA if Fourier *\ \* term turned on, which breaks symmetry. */ if (fparptr->high != NULL) { ptr = fparptr->high; while (ptr != NULL) { if (strncmp (ptr->objtype, "four", 4)==0) mode = 1; ptr = ptr->next; }; }; fparptr->a[10] = pa (mode, fparptr->a[10]); /* PA of 90 and -90 would cause GALFIT to crash prematurely *\ \* because one of the derivatives would be exactly zero. */ if (fparptr->a[10] == 90.) fparptr->a[10] = 89.9999; if (fparptr->a[10] == -90.) fparptr->a[10] = -89.9999; /* Axis ratio can't be too small */ if (!istrunc(fparptr->objtype) && fparptr->a[9] < 0.01) fparptr->a[9] = 0.01; for (i = 1; i<= NPARS; i++) if (fparptr->ia[i] != 1 && fparptr->ia[i] != 0 && fparptr->ia[i] != -1) fparptr->ia[i] = 1; }; fparptr = fparptr->next; }; /************************************************************************\ * Now set the flags for how the type of constraints (strong vs. soft). * * But the constraints are actually linked for the first time in * * check_constraints, which is in mrqmin, because it's easier to deal * * with parameter coupling in the linearized array of fitting params * * instead of in a structure tree, which takes a lot of searching. * \************************************************************************/ while (constr != NULL) { fparptr = fpar; parnum = par2num (constr->parname); if (constr->op == 'o' || constr->op == 'r') { /**********************************************************\ * The first component in the group of constraints is the * * lead, so the ia flag doesn't need to be modified * \**********************************************************/ m = 0; i = 2; while (fparptr != NULL) { m++; if (i <= constr->ncomp && m == constr->comp[i]) { ptr = fparptr; if (strncmp(constr->parname, "f", 1) == 0) { ptr = fparptr->high; while (ptr != NULL && strncmp(ptr->objtype, constr->parname, 1) != 0) ptr = ptr->next; }; i++; if (ptr->ia[parnum] == 0) err = 4; /* Constraining a fixed parameter */ else if (ptr->ia[parnum] > 0 && ptr->ia[parnum] <= 2) ptr->ia[parnum] = 2; else if (ptr->ia[parnum] < 0) ptr->ia[parnum] = -1; }; fparptr = fparptr->next; }; }; constr = constr->next; }; cpar->boxsz[1] = input->convbox[1]; cpar->boxsz[2] = input->convbox[2]; cpar->sampfac = input->sampfac; data->magzpt = input->magzpt; data->dp[1] = input->d[1]; data->dp[2] = input->d[2]; if (data->dp[1] == 0.) data->dp[1] = 1.; if (data->dp[2] == 0.) data->dp[2] = 1.; data->muzpt = data->magzpt + 2.5 * log10(data->dp[1] * data->dp[2]); return (err); }