/***********************************************************************\ * Calculate hybrid transition parameters and constants for derivative * * calculations. * \***********************************************************************/ #include #include #include "structs.h" #include "const.h" void get_hybcons (char *objtype, struct hyblinks *hout, struct hyblinks *hin, struct hybval *hybc1, struct hybval *hybc2, float a[], float kap, float OcosPA, float OsinPA, float *hout_c, float *hin_c, float *hybnorm, float *dfnorm, int *hybshape) { float invp5, pa, rr, rr2, ro_out, ro_in; if (strncmp (objtype, "sersic", 6) == 0) invp5 = 1./a[5]; /********************************************************************\ * Note that hout->type = 1 is where the radial free parameters are * * the break radius and the softening *length*. hout->type = 2 * * is where the free parameters are the break radius and the * * softening *radius*. This is useful because sometimes one may * * want to constrain the outer truncation radius whereas other * * times a user may want to constrain the length. * \********************************************************************/ if (hout->a != NULL) { /* Outer taper */ if (hout->type == 1) { /* Break radius and softening *length* */ hybc1->s = 2.65 - 4.95 * hout->a[3]/hout->a[4]; hybc1->gtr = 2 - hybc1->s; /* has both inner and outer taper. */ if (strncmp (objtype, "sersic", 6) == 0) *dfnorm = -kap*invp5/ a[4] * pow(hout->a[3]/a[4], invp5-1); if (strncmp (objtype, "gauss", 5) == 0) *dfnorm = -hout->a[3] / a[4] / a[4]; if (strncmp (objtype, "expon", 5) == 0) *dfnorm = -1. / a[4]; if (strncmp (objtype, "moffat", 6) == 0) { rr = hout->a[3] / a[4]; *dfnorm = -2 * hout->a[3] / a[4] / a[4] / (1 + rr*rr); }; if (strncmp (objtype, "king", 4) == 0) { rr = hout->a[3] / a[4]; rr2 = a[5]/a[4]; *dfnorm = -2 / rr / a[4] / ( pow (1 + rr*rr, -1./a[6]) - pow (1 + rr2*rr2, -1./a[6]) ) * pow(1+rr*rr, 1./a[6] - 1); }; if (hin->a != NULL && hin->ia[4] == 1) hybc1->dfnorm = *dfnorm; else hybc1->dfnorm = 0.; } else { /* Break radius and softening *radius* */ hybc1->s = 2.65 - (2.65 + 2.3) * hout->a[4]/(hout->a[4] - hout->a[3]); hybc1->gtr = 2 - hybc1->s; if (strncmp (objtype, "sersic", 6) == 0) *dfnorm = -kap*invp5/ a[4] * pow(hout->a[4]/a[4], invp5-1); if (strncmp (objtype, "gauss", 5) == 0) *dfnorm = -hout->a[4] / a[4] / a[4]; if (strncmp (objtype, "expon", 5) == 0) *dfnorm = -1. / a[4]; if (strncmp (objtype, "moffat", 6) == 0) { rr = hout->a[4] / a[4]; *dfnorm = -2 * hout->a[4] / a[4] / a[4] / (1 + rr*rr); }; if (strncmp (objtype, "king", 4) == 0) { rr = hout->a[4] / a[4]; rr2 = a[5]/a[4]; *dfnorm = -2 / rr / a[4] / ( pow (1 + rr*rr, -1./a[6]) - pow (1 + rr2*rr2, -1./a[6]) ) * pow(1+rr*rr, 1./a[6] - 1); }; if (hin->a != NULL && hin->ia[4] == 1) hybc1->dfnorm = *dfnorm; else hybc1->dfnorm = 0.; }; /*********************************************************************\ * If the hybrid doesn't have independent q and PA, then the default * * is to assume that they are the same as the Sersic component * \*********************************************************************/ hout->a[1] = a[1]; /* This is still needed for convolution box */ hout->a[2] = a[2]; /* PA and axis ratio are free parameters */ if (hout->ia[9] != -1) { pa = -hout->a[10] / 180. * PI + PI/2.; hybc1->sinPA = sin (pa); hybc1->cosPA = cos (pa); *hybshape = 1; } else { hout->a[9] = a[9]; hout->a[10] = a[10]; hybc1->sinPA = OsinPA; hybc1->cosPA = OcosPA; }; /************************\ * Diskyness & boxyness * \************************/ if (hout->c0 != NULL) *hout_c = hout->c0->a[1] + 2; else *hout_c = 2.; }; if (hin->a != NULL) { /* Inner taper */ if (hin->type == 1) { /* Break radius and softening *length* */ hybc2->s = 2.65 - 4.95 * hin->a[3]/hin->a[4]; hybc2->gtr = 2 - hybc2->s; if (strncmp (objtype, "sersic", 6) == 0) *dfnorm = -kap*invp5 / a[4] * pow(hin->a[3]/a[4], invp5-1); if (strncmp (objtype, "gauss", 5) == 0) *dfnorm = -hin->a[3] / a[4] / a[4]; if (strncmp (objtype, "expon", 5) == 0) *dfnorm = -1. / a[4]; if (strncmp (objtype, "moffat", 6) == 0) { rr = hin->a[3] / a[4]; *dfnorm = -2 * hin->a[3] / a[4] / a[4] / (1 + rr*rr); }; if (strncmp (objtype, "king", 4) == 0) { rr = hin->a[3] / a[4]; rr2 = a[5]/a[4]; *dfnorm = -2 / rr / a[4] / ( pow (1 + rr*rr, -1./a[6]) - pow (1 + rr2*rr2, -1./a[6]) ) * pow(1+rr*rr, 1./a[6] - 1); }; if (hin->ia[4] == 1) /* If the delta parameter varies... */ hybc2->dfnorm = *dfnorm; else hybc2->dfnorm = 0.; } else { /* Break radius and softening *radius* */ hybc2->s = 2.65 - (2.65 + 2.3) * hin->a[4]/(hin->a[4] - hin->a[3]); hybc2->gtr = 2 - hybc2->s; if (strncmp (objtype, "sersic", 6) == 0) *dfnorm = -kap*invp5 / a[4] * pow(hin->a[4]/a[4], invp5-1); if (strncmp (objtype, "gauss", 5) == 0) *dfnorm = -hin->a[4] / a[4] / a[4]; if (strncmp (objtype, "expon", 5) == 0) *dfnorm = -1. / a[4]; if (strncmp (objtype, "moffat", 6) == 0) { rr = hin->a[4] / a[4]; *dfnorm = -2 * hin->a[4] / a[4] / a[4] / (1 + rr*rr); }; if (strncmp (objtype, "king", 4) == 0) { rr = hin->a[4] / a[4]; rr2 = a[5]/a[4]; *dfnorm = -2 / rr / a[4] / ( pow (1 + rr*rr, -1./a[6]) - pow (1 + rr2*rr2, -1./a[6]) ) * pow(1+rr*rr, 1./a[6] - 1); }; if (hin->ia[4] == 1) hybc2->dfnorm = *dfnorm; else hybc2->dfnorm = 0.; }; /*********************************************************************\ * If the hybrid doesn't have independent q and PA, then the default * * is to assume that they are the same as the Sersic component. * \*********************************************************************/ hin->a[1] = a[1]; /* This is still needed for convolution box */ hin->a[2] = a[2]; /* PA and axis ratio are free parameters */ if (hin->ia[9] != -1) { pa = -hin->a[10] / 180. * PI + PI/2.; hybc2->sinPA = sin (pa); hybc2->cosPA = cos (pa); *hybshape = 1; } else { hin->a[9] = a[9]; hin->a[10] = a[10]; hybc2->sinPA = OsinPA; hybc2->cosPA = OcosPA; }; /*****************************************************************\ * If there's an inner taper (hin->type=1), then the flux is * * normalized at the outer break radius. If there's an outer * * taper (hin->type=2), then the flux is normalized at the inner * * break radius * \*****************************************************************/ if (hin->type == 1) { if (strncmp (objtype, "sersic", 6) == 0) { rr = pow(hin->a[3]/a[4], invp5); *hybnorm = exp(-kap*(rr-1.)); }; if (strncmp (objtype, "gauss", 5) == 0) { rr = pow(hin->a[3]/a[4], 2); *hybnorm = exp(-rr / 2.); }; if (strncmp (objtype, "expon", 5) == 0) *hybnorm = exp(- hin->a[3]/a[4]); if (strncmp (objtype, "moffat", 6) == 0) { rr = hin->a[4]/a[4]; *hybnorm = 1./ pow(1 + rr*rr, a[5]); }; } else { if (strncmp (objtype, "sersic", 6) == 0) { rr = pow(hin->a[4]/a[4], invp5); *hybnorm = exp(-kap*(rr-1.)); }; if (strncmp (objtype, "gauss", 5) == 0) { rr = pow(hin->a[4]/a[4], 2); *hybnorm = exp(-rr / 2.); }; if (strncmp (objtype, "expon", 5) == 0) *hybnorm = exp(- hin->a[4]/a[4]); if (strncmp (objtype, "moffat", 6) == 0) { rr = hin->a[4]/a[4]; *hybnorm = 1./ pow(1 + rr*rr, a[5]); }; }; /************************\ * Diskyness & boxyness * \************************/ if (hin->c0 != NULL) *hin_c = hin->c0->a[1] + 2; else *hin_c = 2.; }; }