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