Changeset 8882 for trunk/psphot/src/models/pmModel_RGAUSS.c
- Timestamp:
- Sep 22, 2006, 2:29:31 AM (20 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/models/pmModel_RGAUSS.c (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/models/pmModel_RGAUSS.c
r4954 r8882 20 20 psF32 *PAR = params->data.F32; 21 21 22 psF32 X = x->data.F32[0] - PAR[ 2];23 psF32 Y = x->data.F32[1] - PAR[ 3];24 psF32 px = PAR[ 4]*X;25 psF32 py = PAR[ 5]*Y;26 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[ 6]*X*Y;22 psF32 X = x->data.F32[0] - PAR[PM_PAR_XPOS]; 23 psF32 Y = x->data.F32[1] - PAR[PM_PAR_YPOS]; 24 psF32 px = PAR[PM_PAR_SXX]*X; 25 psF32 py = PAR[PM_PAR_SYY]*Y; 26 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; 27 27 28 28 // psF32 FACT = 1 + 5*exp(-5*PAR[7]); 29 29 30 psF32 p = pow(z, PAR[ 7] - 1.0);30 psF32 p = pow(z, PAR[PM_PAR_7] - 1.0); 31 31 psF32 r = 1.0 / (1 + z + z*p); 32 psF32 f = PAR[ 1]*r + PAR[0];32 psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; 33 33 34 34 if (deriv != NULL) { 35 35 // note difference from a pure gaussian: q = params->data.F32[1]*r 36 psF32 t = PAR[ 1]*r*r;37 psF32 q = t*(1 + PAR[ 7]*p);36 psF32 t = PAR[PM_PAR_FLUX]*r*r; 37 psF32 q = t*(1 + PAR[PM_PAR_7]*p); 38 38 39 39 deriv->data.F32[0] = +1.0; 40 40 deriv->data.F32[1] = +r; 41 deriv->data.F32[2] = q*(2.0*px*PAR[ 4] + PAR[6]*Y);42 deriv->data.F32[3] = q*(2.0*py*PAR[ 5] + PAR[6]*X);41 deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); 42 deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); 43 43 deriv->data.F32[4] = -q*px*X*2; 44 44 deriv->data.F32[5] = -q*py*Y*2; … … 60 60 psF32 *PAR = params->data.F32; 61 61 62 psF64 A1 = PS_SQR(PAR[ 4]);63 psF64 A2 = PS_SQR(PAR[ 5]);64 psF64 A3 = PS_SQR(PAR[ 6]);62 psF64 A1 = PS_SQR(PAR[PM_PAR_SXX]); 63 psF64 A2 = PS_SQR(PAR[PM_PAR_SYY]); 64 psF64 A3 = PS_SQR(PAR[PM_PAR_SXY]); 65 65 psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); 66 66 // Area is equivalent to 2 pi sigma^2 … … 69 69 norm = 0.0; 70 70 for (z = 0.005; z < 50; z += 0.01) { 71 f = 1.0 / (1 + z + pow(z, PAR[ 7]));71 f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); 72 72 norm += f; 73 73 } … … 87 87 88 88 if (flux <= 0) return (1.0); 89 if (PAR[ 1] <= 0) return (1.0);90 if (flux >= PAR[ 1]) return (1.0);89 if (PAR[PM_PAR_FLUX] <= 0) return (1.0); 90 if (flux >= PAR[PM_PAR_FLUX]) return (1.0); 91 91 92 92 // if Sx == Sy, sigma = Sx == Sy 93 psF64 sigma = hypot (1.0 / PAR[ 4], 1.0 / PAR[5]) / sqrt(2.0);93 psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0); 94 94 psF64 dz = 1.0 / (2.0 * sigma*sigma); 95 psF64 limit = flux / PAR[ 1];95 psF64 limit = flux / PAR[PM_PAR_FLUX]; 96 96 97 97 // we can do this much better with intelligent choices here 98 98 for (z = 0.0; z < 20.0; z += dz) { 99 p = pow(z, PAR[ 7]);99 p = pow(z, PAR[PM_PAR_7]); 100 100 f = 1.0 / (1 + z + p); 101 101 if (f < limit) break;
Note:
See TracChangeset
for help on using the changeset viewer.
