Changeset 8882 for trunk/psphot/src/models/pmModel_QGAUSS.c
- Timestamp:
- Sep 22, 2006, 2:29:31 AM (20 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/models/pmModel_QGAUSS.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/models/pmModel_QGAUSS.c
r5593 r8882 26 26 psF32 *PAR = params->data.F32; 27 27 28 psF32 X = x->data.F32[0] - PAR[ 2];29 psF32 Y = x->data.F32[1] - PAR[ 3];30 psF32 px = PAR[ 4]*X;31 psF32 py = PAR[ 5]*Y;32 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[ 6]*X*Y;28 psF32 X = x->data.F32[0] - PAR[PM_PAR_XPOS]; 29 psF32 Y = x->data.F32[1] - PAR[PM_PAR_YPOS]; 30 psF32 px = PAR[PM_PAR_SXX]*X; 31 psF32 py = PAR[PM_PAR_SYY]*Y; 32 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; 33 33 34 psF32 r = 1.0 / (1 + PAR[ 7]*z + pow(z, QG_S1));35 psF32 f = PAR[ 1]*r + PAR[0];34 psF32 r = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1)); 35 psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; 36 36 37 37 if (deriv != NULL) { 38 38 // note difference from a pure gaussian: q = params->data.F32[1]*r 39 psF32 t = PAR[ 1]*r*r;40 psF32 q = t*(PAR[ 7] + QG_S1*pow(z, dQG_S1));39 psF32 t = PAR[PM_PAR_FLUX]*r*r; 40 psF32 q = t*(PAR[PM_PAR_7] + QG_S1*pow(z, dQG_S1)); 41 41 42 42 deriv->data.F32[0] = +1.0; 43 43 deriv->data.F32[1] = +r; 44 deriv->data.F32[2] = q*(2.0*px*PAR[ 4] + PAR[6]*Y);45 deriv->data.F32[3] = q*(2.0*py*PAR[ 5] + PAR[6]*X);44 deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); 45 deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); 46 46 deriv->data.F32[4] = -2.0*q*px*X; 47 47 deriv->data.F32[5] = -2.0*q*py*Y; … … 111 111 psF32 *PAR = params->data.F32; 112 112 113 psF64 A1 = PS_SQR(PAR[ 4]);114 psF64 A2 = PS_SQR(PAR[ 5]);115 psF64 A3 = PS_SQR(PAR[ 6]);113 psF64 A1 = PS_SQR(PAR[PM_PAR_SXX]); 114 psF64 A2 = PS_SQR(PAR[PM_PAR_SYY]); 115 psF64 A3 = PS_SQR(PAR[PM_PAR_SXY]); 116 116 psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); 117 117 // Area is equivalent to 2 pi sigma^2 … … 120 120 norm = 0.0; 121 121 for (z = 0.005; z < 50; z += 0.01) { 122 f = 1.0 / (1 + PAR[ 7]*z + pow(z, QG_S1));122 f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1)); 123 123 norm += f; 124 124 } 125 125 norm *= 0.01; 126 126 127 psF64 Flux = PAR[ 1] * Area * norm;127 psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm; 128 128 129 129 return(Flux); … … 138 138 139 139 if (flux <= 0) return (1.0); 140 if (PAR[ 1] <= 0) return (1.0);141 if (flux >= PAR[ 1]) return (1.0);140 if (PAR[PM_PAR_FLUX] <= 0) return (1.0); 141 if (flux >= PAR[PM_PAR_FLUX]) return (1.0); 142 142 143 143 // if Sx == Sy, sigma = Sx == Sy 144 psF64 sigma = hypot (1.0 / PAR[ 4], 1.0 / PAR[5]) / sqrt(2.0);144 psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0); 145 145 psF64 dz = 1.0 / (2.0 * sigma*sigma); 146 psF64 limit = flux / PAR[ 1];146 psF64 limit = flux / PAR[PM_PAR_FLUX]; 147 147 148 148 // we can do this much better with intelligent choices here 149 149 for (z = 0.0; z < 20.0; z += dz) { 150 f = 1.0 / (1 + PAR[ 7]*z + pow(z, QG_S1));150 f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1)); 151 151 if (f < limit) break; 152 152 } … … 184 184 185 185 dP = 0; 186 dP += PS_SQR(dPAR[ 4] / PAR[4]);187 dP += PS_SQR(dPAR[ 5] / PAR[5]);186 dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); 187 dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); 188 188 dP = sqrt (dP); 189 189 190 190 status = true; 191 191 status &= (dP < 0.5); 192 status &= (PAR[ 1] > 0);193 status &= ((dPAR[ 1]/PAR[1]) < 0.5);192 status &= (PAR[PM_PAR_FLUX] > 0); 193 status &= ((dPAR[PM_PAR_FLUX]/PAR[PM_PAR_FLUX]) < 0.5); 194 194 195 195 if (status) return true;
Note:
See TracChangeset
for help on using the changeset viewer.
