Changeset 8882 for trunk/psphot/src/models/pmModel_ZGAUSS.c
- Timestamp:
- Sep 22, 2006, 2:29:31 AM (20 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/models/pmModel_ZGAUSS.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/models/pmModel_ZGAUSS.c
r4954 r8882 23 23 psF32 *PAR = params->data.F32; 24 24 25 psF32 X = x->data.F32[0] - PAR[ 2];26 psF32 Y = x->data.F32[1] - PAR[ 3];27 psF32 px = PAR[ 4]*X;28 psF32 py = PAR[ 5]*Y;29 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[ 6]*X*Y;25 psF32 X = x->data.F32[0] - PAR[PM_PAR_XPOS]; 26 psF32 Y = x->data.F32[1] - PAR[PM_PAR_YPOS]; 27 psF32 px = PAR[PM_PAR_SXX]*X; 28 psF32 py = PAR[PM_PAR_SYY]*Y; 29 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; 30 30 31 31 psF32 pr = PAR8*z; 32 psF32 p = pow(z, PAR[ 7] - 1.0);32 psF32 p = pow(z, PAR[PM_PAR_7] - 1.0); 33 33 psF32 r = 1.0 / (1 + z*p + SQ(SQ(pr))); 34 psF32 f = PAR[ 1]*r + PAR[0];34 psF32 f = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY]; 35 35 36 36 if (deriv != NULL) { 37 37 // note difference from a pure gaussian: q = params->data.F32[1]*r 38 psF32 t = PAR[ 1]*r*r;39 psF32 q = t*(PAR[ 7]*p + 4*PAR8*pr*pr*pr);38 psF32 t = PAR[PM_PAR_FLUX]*r*r; 39 psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR8*pr*pr*pr); 40 40 41 41 deriv->data.F32[0] = +1.0; 42 42 deriv->data.F32[1] = +r; 43 deriv->data.F32[2] = q*(2.0*px*PAR[ 4] + PAR[6]*Y);44 deriv->data.F32[3] = q*(2.0*py*PAR[ 5] + PAR[6]*X);43 deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y); 44 deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X); 45 45 deriv->data.F32[4] = -q*px*X; 46 46 deriv->data.F32[5] = -q*py*Y; … … 57 57 psF32 *PAR = params->data.F32; 58 58 59 psF64 A1 = PS_SQR(PAR[ 4]);60 psF64 A2 = PS_SQR(PAR[ 5]);61 psF64 A3 = PS_SQR(PAR[ 6]);59 psF64 A1 = PS_SQR(PAR[PM_PAR_SXX]); 60 psF64 A2 = PS_SQR(PAR[PM_PAR_SYY]); 61 psF64 A3 = PS_SQR(PAR[PM_PAR_SXY]); 62 62 psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); 63 63 // Area is equivalent to 2 pi sigma^2 … … 67 67 psF32 pr = PAR8*z; 68 68 for (z = 0.005; z < 50; z += 0.01) { 69 f = 1.0 / (1 + pow(z, PAR[ 7]) + SQ(SQ(pr)));69 f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr))); 70 70 norm += f; 71 71 } 72 72 norm *= 0.01; 73 73 74 psF64 Flux = PAR[ 1] * Area * norm;74 psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm; 75 75 76 76 return(Flux); … … 89 89 90 90 if (flux <= 0) return (1.0); 91 if (PAR[ 1] <= 0) return (1.0);92 if (flux >= PAR[ 1]) return (1.0);91 if (PAR[PM_PAR_FLUX] <= 0) return (1.0); 92 if (flux >= PAR[PM_PAR_FLUX]) return (1.0); 93 93 94 94 // convert Sx,Sy,Sxy to major/minor axes 95 shape.sx = 1.0 / PAR[ 4];96 shape.sy = 1.0 / PAR[ 5];97 shape.sxy = PAR[ 6];95 shape.sx = 1.0 / PAR[PM_PAR_SXX]; 96 shape.sy = 1.0 / PAR[PM_PAR_SYY]; 97 shape.sxy = PAR[PM_PAR_SXY]; 98 98 99 99 axes = EllipseShapeToAxes (shape); 100 100 psF64 dr = 1.0 / axes.major; 101 psF64 limit = flux / PAR[ 1];101 psF64 limit = flux / PAR[PM_PAR_FLUX]; 102 102 103 103 // XXX : we can do this faster with an intelligent starting choice … … 105 105 z = SQ(r); 106 106 pr = PAR8*z; 107 f = 1.0 / (1 + pow(z, PAR[ 7]) + SQ(SQ(pr)));107 f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr))); 108 108 if (f < limit) break; 109 109 }
Note:
See TracChangeset
for help on using the changeset viewer.
