/***************************************************************************\ * Calculate truncation transition parameters and constants for derivative * * calculations for length and height truncations, as opposed to radial * * truncations that are calculated in "get_trunc_cons.c" * \***************************************************************************/ #include #include #include "structs.h" #include "nrutil.h" #include "const.h" #define EPS 1e-5 float bessk1 (float x); float bessk0 (float x); void get_trunc_cons2 (int ttype, struct trunc_links *ht_out, struct trunc_links *ht_in, struct trunc_val *ht_out_cons, struct trunc_val *ht_in_cons, struct trunc_links *lg_out, struct trunc_links *lg_in, struct trunc_val *lg_out_cons, struct trunc_val *lg_in_cons, float a[], float *trunc_norm, float *ht_dfnorm, float *lg_dfnorm) { float hh, rr, K0, K1, sech2, ro_out, ro_in; /**********************************************************************\ * Note that trunc_out->type = 1 is where the radial free parameters * * are the break radius and the softening *length*. * * trunc_out->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. * * * * Note also that outer taper is normalized at the center, except * * when the model has both an inner and outer truncation. In that * * case, the flux is normalized at the inner truncation radius. * \**********************************************************************/ if (ttype == 100 || ttype == 111 || ttype == 101 || ttype == 110) *trunc_norm = 1.; if (ht_out != NULL && ht_out->a != NULL) { /* Height outer truncation */ /*--------------------------------------------------------------\ | For outer truncation, the flux is normalized at the center | | unless there's also an inner truncation. In the latter | | scenario, the flux is normalized at the break radius | | calculated in the section far below on "Inner taper". | \--------------------------------------------------------------*/ if (ht_out->lrtype == 1) { /* Break radius and softening *length* */ ht_out_cons->s = 2.65 - 4.95 * ht_out->a[4]/ht_out->a[5]; ht_out_cons->gtr = 2 - ht_out_cons->s; ht_out_cons->dfnorm = 0.; } else { /* Break radius and softening *radius* */ ht_out_cons->s = 2.65 - (2.65 + 2.3) * ht_out->a[5]/ (ht_out->a[5] - ht_out->a[4]); ht_out_cons->gtr = 2 - ht_out_cons->s; ht_out_cons->dfnorm = 0.; }; /*--------------------------------------------------------------------\ | If the truncation doesn't have independent q and PA, then the | | default is to assume that they are the same as the Sersic comp.. | \--------------------------------------------------------------------*/ if (ht_out->ia[1] == -1) { ht_out->a[1] = a[1]; /* This is still needed for convolution box */ ht_out->a[2] = a[2]; }; }; if (lg_out != NULL && lg_out->a != NULL) { /* Length outer truncation */ if (lg_out->lrtype == 1) { /* Break radius and softening *length* */ lg_out_cons->s = 2.65 - 4.95 * lg_out->a[4]/lg_out->a[5]; lg_out_cons->gtr = 2 - lg_out_cons->s; lg_out_cons->dfnorm = 0.; } else { /* Break radius and softening *radius* */ lg_out_cons->s = 2.65 - (2.65 + 2.3) * lg_out->a[5]/ (lg_out->a[5] - lg_out->a[4]); lg_out_cons->gtr = 2 - lg_out_cons->s; lg_out_cons->dfnorm = 0.; }; if (lg_out->ia[1] == -1) { lg_out->a[1] = a[1]; /* This is still needed for convolution box */ lg_out->a[2] = a[2]; }; }; /**********************************************************************\ * INNER TRUNCATION on a function. If there is an inner truncation * * then the flux is normalized to the break radius, near which the * * product of the hyperbolic tangent and the profile function reaches * * a maximum. * \**********************************************************************/ if (ttype != 100 && ttype != 101 && ttype != 110 && ttype != 111 ) { /*----------------------------------------------------------------\ | Now the truncation normalization and derivatives with respect | | to the break radius, "dfnorm." | \*---------------------------------------------------------------*/ if (ht_in != NULL) { /* Height inner truncation */ if (ht_in->lrtype == 1) { /* Break radius and softening *length* */ ht_in_cons->s = 2.65 - 4.95 * ht_in->a[4] / ht_in->a[5]; hh = ht_in->a[4] / a[4]; } else { /* Break radius and softening *radius* */ ht_in_cons->s = 2.65 - (2.65 + 2.3) * ht_in->a[5] / (ht_in->a[5] - ht_in->a[4]); hh = ht_in->a[5] / a[4]; }; ht_in_cons->gtr = 2 - ht_in_cons->s; sech2 = 1./pow(cosh(hh), 2.); *ht_dfnorm = -sinh(hh)/cosh(hh) * 2 / a[4]; }; if (lg_in != NULL) { /* Length inner truncation */ if (lg_in->lrtype == 1) { /* Break radius and softening *length* */ lg_in_cons->s = 2.65 - 4.95 * lg_in->a[4] / lg_in->a[5]; rr = lg_in->a[4] / a[5]; } else { /* Break radius and softening *radius* */ lg_in_cons->s = 2.65 - (2.65 + 2.3) * lg_in->a[5] / (lg_in->a[5] - lg_in->a[4]); rr = lg_in->a[5] / a[5]; }; lg_in_cons->gtr = 2 - lg_in_cons->s; K1 = bessk1 (rr); K0 = bessk0 (rr); *lg_dfnorm = (K1 + rr * (-K0 - K1/rr))/ K1/rr/a[5]; }; if (ttype == 120 || ttype == 121 || ttype == 130 || ttype == 131) *trunc_norm = sech2; else if (ttype == 102 || ttype == 103 || ttype == 112 || ttype == 113) *trunc_norm = rr * K1; else if (ttype == 122 || ttype == 123 || ttype == 132 || ttype == 133) *trunc_norm = rr * K1 * sech2; if (ht_in != NULL) { if ((ht_in->lrtype == 1 && ht_in->ia[3] == 1) || (ht_in->lrtype == 2 && ht_in->ia[4] == 1)) ht_in_cons->dfnorm = *ht_dfnorm; else ht_in_cons->dfnorm = 0.; } else ht_in_cons->dfnorm = 0.; if (lg_in != NULL) { if ((lg_in->lrtype == 1 && lg_in->ia[3] == 1) || (lg_in->lrtype == 2 && lg_in->ia[4] == 1)) lg_in_cons->dfnorm = *lg_dfnorm; else lg_in_cons->dfnorm = 0.; } else lg_in_cons->dfnorm = 0.; if (ht_in != NULL && ht_in->ia[1] == -1) { ht_in->a[1] = a[1]; /* This is still needed for convolution box */ ht_in->a[2] = a[2]; }; if (lg_in != NULL && lg_in->ia[1] == -1) { lg_in->a[1] = a[1]; /* This is still needed for convolution box */ lg_in->a[2] = a[2]; }; }; }