/***************************************************************************\ * This subroutine is controls the sanity of the parameter step sizes. * * Sometimes parameter step size can be way too big even though the * * sentiment is in the right direction. This subroutine enforces changes * * that are not too big. * \***************************************************************************/ #include #include #include "nrutil.h" #include "structs.h" #define DELMAG 10. /* max change in mag */ #define DELAX 100. /* max change in axis ratio */ #define MAXLENGRATIO 10. /* max change in scale length ratio */ #define MAXNRATIO 10. /* maximum change in n (ratio) */ #define MAXDPA 500. /* maximum change in PA (degrees) */ #define MAXDELTA 100. /* maximum change in position (pix) */ #define MINAXRATIO 0.01 /* minimum allwed axis ratio */ #define MINLENG 0.01 /* minimum allowed scale length */ float pa (int mode, float instrpa); struct findpar { int mode; int out; }; int istrunc (char *objtype); void par_monitor (int ma, float a[], int ia[], float atry[], double da[], struct fitpars *fpar) { int matchpar (int l, struct fitpars *fpar, struct findpar *k, char obj[]); extern int (*pfunc)(const char *, ...); extern char device[]; int j, l; struct findpar k = {0, 0}; char obj[10]; for (j=0, l=1; l<=ma; l++) { if (ia[l] == 1) { ++j; atry[l] = a[l] + da[j]; /*--------------------------------------------------------\ | Find the parameter corresponding to the l-th element, | | a[l]. Return the matched parameter in k. | \--------------------------------------------------------*/ matchpar (l, fpar, &k, obj); /*------------------\ | Position angle | \------------------*/ if (k.out == 10) atry[l] = pa (k.mode, atry[l]); /*----------------------------------------------\ | AXIS RATIO: Constrain to be between 0 and 1 | | have to be careful what to do here, because | | it affects convergence of bending, fourier, | | parameters. | \----------------------------------------------*/ if (k.out == 9) { if (atry[l] > 1.) { if (1./atry[l] < a[l]) atry[l] = 0.99; else /* Invert ax-ratio */ atry[l] = 1./atry[l]; } if (atry[l] < 0.) /* Make sure the axis ratio */ atry[l] = 0.05; /* is not less than 0. */ }; /*----------------------------------------------------\ * SCALELENGTHS: Constrain to be greater than 1e-2 | \----------------------------------------------------*/ if (k.out == 4) { /* Make sure length doesn't change so much at once */ /* if (atry[l] / a[l] >= MAXLENGRATIO) atry[l] = a[l] * MAXLENGRATIO; else if (a[l] / atry[l] >= MAXLENGRATIO) atry[l] = a[l] / MAXLENGRATIO; */ if ( a[l] <= 0.00 ) a[l] = MINLENG; if ( atry[l] <= 0.00 ) atry[l] = MINLENG; }; /*-----------------------------------------------------\ | POWERLAW: Constrain alpha and n to be less than 20 | | and greater than 0.01. | \-----------------------------------------------------*/ if (k.out==5 && strncmp (obj, "edgedisk", 8) != 0) { if (strncmp (obj, "king", 4) != 0) { /* if (fabs(atry[l] / a[l]) >= MAXNRATIO) atry[l] = a[l] * MAXNRATIO; else if (a[l] / atry[l] >= MAXNRATIO) atry[l] = a[l] / MAXNRATIO; */ if (atry[l] > 20. && !istrunc(obj)) atry[l] = 1e30; else if (atry[l] < 0.0) atry[l] = 1e30; }; }; /*-----------------------------------------------------------\ | PA, and (X, Y) POSITIONS: can be allowed to be negative, | | but if the change is too big, the initial guess is | | probably way too bad. Hold it fixed for now.... | \-----------------------------------------------------------*/ if (k.out==1 || k.out==2) { if (atry[l] - a[l] >= MAXDELTA) atry[l] = a[l] + MAXDELTA; else if (a[l] - atry[l] >= MAXDELTA) atry[l] = a[l] - MAXDELTA; }; /*-------------------------------------------------------------\ | DISKYNESS/BOXYNESS: constrain to between -2 and positive 2 | \-------------------------------------------------------------*/ /* if (l % 10 == 0) { if (atry[l] <= -2.0) atry[l] = -2.0; }; */ /*-------------------------------------------\ | Constrain everything else to be positive | \-------------------------------------------*/ if (k.out >= 6 && k.out <= 8) if (atry[l] < 0.) atry[l] = 0.01; }; }; } /****************************************************************************/ int matchpar (int l, struct fitpars *fpar, struct findpar *kout, char obj[]) { int i, mode, found=0; struct fitpars *ptr; static int k=0, high = 0; while (fpar != NULL && !found ) { /*------------------------------------------------------------------\ | Watch out for truncation parameters. Those parameters that are | | in the truncation structures count, but not those in the *high* | | structures, which are only pointers to the truncation structs. | \------------------------------------------------------------------*/ if (!istrunc (fpar->objtype) || (istrunc (fpar->objtype) && high == 0)) { for (i=1; i <= fpar->npars; i++) { if (!found) { k++; if (k==l) { found = 1; strncpy (obj, fpar->objtype, 10); if (high == 0 && strncmp (fpar->objtype, "sky", 3) != 0) { mode = 2; k = i; /*-----------------------------------------------\ | If Fourier term is turned on, the classic PA | | ranges from 0 to 360 instead of -90 to 90, | | since the Fourier modes break symmetry. | \-----------------------------------------------*/ if (i == 10 && fpar->high != NULL) { ptr = fpar->high; while (ptr != NULL) { if (strncmp (ptr->objtype, "four", 4)==0) mode = 1; ptr = ptr->next; }; }; } else if (strncmp (fpar->objtype, "sky", 3) == 0) k = 0; if (high == 1) { /* Check phase angle for fourier term */ if (strncmp (fpar->objtype, "four", 4)==0 && i%2 == 0) { k = 10; mode = i/2; /* Check the inclination angle and sky PA */ /* for coord rot. */ } else if ((strncmp (fpar->objtype, "power", 5)==0 || strncmp (fpar->objtype, "log", 3)==0 ) && (i==9 || i==10)) { k = 10; mode = 1; }; }; }; }; }; }; /* Go down high order list if high == 1 */ if (fpar->high != NULL && !found && high == 0) { high += 1; found = matchpar (l, fpar->high, kout, obj); high -= 1; /* set high = 0 once out of the high-order list */ }; fpar = fpar->next; }; if (found == 1) { kout->out = k; kout->mode = mode; k = 0; /* Need to set k back to 0 because it's static. */ found = 2; /* And set found to 2 so subsequent returns back *\ * to the recursive loop doesn't enter in here * \* again. */ } else if (found == 0 && high == 0 && fpar == NULL ) k = 0; return (found); }