00001
00002
00003
00004
00005
00006
00007 #include <stdio.h>
00008 #include <math.h>
00009
00010 #include "imcore.h"
00011 #include "imcore_radii.h"
00012 #include "floatmath.h"
00013 #include "util.h"
00014 #include "ap.h"
00015
00016 static float fraction (float x, float y, float r_out);
00017 static void dchole (double a[IMNUM+1][IMNUM+1], double b[IMNUM+1], int n);
00018
00019 extern float imcore_exprad(float thresh, float peak, float areal0,
00020 float rcores[], int naper) {
00021 float pk,r_t,rad;
00022
00023
00024
00025 pk = MAX(1.5*thresh,peak);
00026 r_t = sqrtf(areal0/M_PI);
00027 rad = 5.0*r_t/logf(pk/thresh);
00028 rad = MAX(r_t,MIN(5.0*r_t,MIN(rad,rcores[naper-1])));
00029 return(rad);
00030 }
00031
00032 extern float imcore_kronrad(float areal0, float rcores[], float cflux[],
00033 int naper) {
00034 int i;
00035 float r_t,rad,wt;
00036
00037
00038
00039 r_t = sqrtf(areal0/M_PI);
00040 rad = 0.5*rcores[0]*cflux[0];
00041 for (i = 1; i < naper; i++) {
00042 wt = MAX(0.0,cflux[i]-cflux[i-1]);
00043 rad += 0.5*(rcores[i] + rcores[i-1])*wt;
00044 }
00045 rad /= cflux[MIN(naper-1,7)];
00046 rad = MAX(r_t,MIN(5.0*r_t,MIN(2.0*rad,rcores[naper-1])));
00047 return(rad);
00048 }
00049
00050 extern float imcore_petrad (float areal0, float rcores[], float cflux[],
00051 int naper) {
00052 int j;
00053 float eta,r_t,etaold,r1,r2,r3,r4,r5,r_petr;
00054
00055
00056
00057 r_t = sqrtf(areal0/M_PI);
00058 eta = 1.0;
00059 etaold = eta;
00060 j = 1;
00061 while (eta > 0.2 && j < naper) {
00062 etaold = eta;
00063 r1 = rcores[j]*rcores[j]/(rcores[j-1]*rcores[j-1]) - 1.0;
00064 r2 = cflux[j]/cflux[j-1] - 1.0;
00065 eta = r2/r1;
00066 j++;
00067 }
00068 if (j == naper) {
00069 r_petr = rcores[naper-1];
00070 } else {
00071 r1 = rcores[j]*rcores[j];
00072 r2 = rcores[j-1]*rcores[j-1];
00073 r3 = rcores[j-2]*rcores[j-2];
00074 r4 = (etaold - 0.2)/(etaold - eta);
00075 r5 = (0.2 - eta)/(etaold - eta);
00076 r_petr = r4*sqrt(0.5*(r1 + r2)) + r5*sqrt(0.5*(r2 + r3));
00077 }
00078 r_petr = MAX(r_t,MIN(5.0*r_t,MIN(2.0*r_petr,rcores[naper-1])));
00079 return(r_petr);
00080 }
00081
00082 extern void imcore_flux(ap_t *ap, float parm[IMNUM][NPAR], int nbit,
00083 float apertures[], float cflux[]) {
00084 double aa[IMNUM+1][IMNUM+1],bb[IMNUM+1];
00085 float *map,parrad[IMNUM],cn[IMNUM],xi,yi,d,a2j,a2i,di,dj,argi,argj;
00086 float xmin,xmax,ymin,ymax,t,xj,yj,tj,xk,yk,tk;
00087 unsigned char *mflag;
00088 long nx,ny;
00089 int i,ix1,ix2,iy1,iy2,ii,kk,j,k;
00090
00091
00092
00093 map = ap->indata;
00094 mflag = ap->mflag;
00095 nx = ap->lsiz;
00096 ny = ap->csiz;
00097
00098
00099
00100 for (i = 0; i < nbit; i++) {
00101 parrad[i] = apertures[i] + 0.5;
00102 cn[i] = 1.0/(M_PI*apertures[i]*apertures[i]);
00103 }
00104
00105
00106
00107 for(i = 0; i < nbit; i++) {
00108 aa[i][i] = cn[i];
00109 if (nbit > 1) {
00110 xi = parm[i][1];
00111 yi = parm[i][2];
00112 for(j = i+1; j < nbit; j++) {
00113 d = sqrtf((xi-parm[j][1])*(xi-parm[j][1])
00114 + (yi-parm[j][2])*(yi-parm[j][2]));
00115 if(d >= apertures[i]+apertures[j]) {
00116 aa[j][i] = 0.0;
00117 aa[i][j] = aa[j][i];
00118 } else {
00119 a2i = apertures[i]*apertures[i];
00120 a2j = apertures[j]*apertures[j];
00121 di = 0.5*(d + (a2i - a2j)/d);
00122 dj = 0.5*(d - (a2i - a2j)/d);
00123 argi = di/apertures[i];
00124 argj = dj/apertures[j];
00125 if (di < 0.0) {
00126 argi = 1.0;
00127 argj = 0.0;
00128 }
00129 if (dj < 0.0) {
00130 argi = 0.0;
00131 argj = 1.0;
00132 }
00133 aa[j][i] = cn[i]*cn[j]*
00134 (a2i*(acosf(argi) - argi*(sqrtf(1.0 - argi*argi))) +
00135 a2j*(acosf(argj) - argj*(sqrtf(1.0 - argj*argj))));
00136 aa[i][j] = aa[j][i];
00137 }
00138 }
00139 }
00140 }
00141
00142
00143
00144 for (i = 0; i < nbit; i++)
00145 bb[i] = 0.0;
00146
00147
00148
00149 xmin = 1.0e6;
00150 xmax = -1.0e6;
00151 ymin = 1.0e6;
00152 ymax = -1.0e6;
00153 for(i = 0; i < nbit; i++) {
00154 xi = parm[i][1];
00155 yi = parm[i][2];
00156 xmin = MIN(xmin, xi - parrad[i]);
00157 xmax = MAX(xmax, xi + parrad[i]);
00158 ymin = MIN(ymin, yi - parrad[i]);
00159 ymax = MAX(ymax, yi + parrad[i]);
00160 }
00161 ix1 = MAX(0,(int)xmin-1);
00162 ix2 = MIN(nx-1,(int)xmax);
00163 iy1 = MAX(0,(int)ymin-1);
00164 iy2 = MIN(ny-1,(int)ymax);
00165
00166
00167
00168
00169 for(ii = iy1; ii <= iy2; ii++) {
00170 kk = ii*nx;
00171 for(i = ix1; i <= ix2; i++) {
00172 if (mflag[kk+i] == MF_CLEANPIX || mflag[kk+i] == MF_OBJPIX) {
00173 t = map[kk+i];
00174 for(j = 0; j < nbit; j++) {
00175 xj = i - parm[j][1] + 1.0;
00176 yj = ii - parm[j][2] + 1.0;
00177 bb[j] += fraction(xj,yj,apertures[j])*t;
00178 }
00179 } else {
00180 for (j = 0; j < nbit; j++) {
00181 xj = i - parm[j][1] + 1.0;
00182 yj = ii - parm[j][2] + 1.0;
00183 tj = fraction(xj,yj,apertures[j]);
00184 aa[j][j] -= tj*tj*cn[j]*cn[j];
00185 for (k = j + 1; k < nbit; k++) {
00186 xk = i - parm[k][1] + 1.0;
00187 yk = ii - parm[k][2] + 1.0;
00188 tk = fraction(xk,yk,apertures[k]);
00189 aa[j][k] -= tk*tj*cn[j]*cn[k];
00190 aa[k][j] = aa[j][k];
00191 }
00192 }
00193 }
00194 }
00195 }
00196
00197
00198
00199 if (nbit > 1) {
00200 for (k = 0; k < nbit; k++) {
00201 for (j = k + 1; j < nbit; j++) {
00202 aa[j][k] = 0.99*MIN(aa[k][k],aa[j][j]);
00203 aa[k][j] = aa[j][k];
00204 }
00205 }
00206 }
00207
00208
00209
00210 if (nbit == 1) {
00211 cflux[0] = bb[0];
00212 return;
00213 }
00214
00215
00216
00217 dchole(aa,bb,nbit);
00218 for(i = 0; i < nbit; i++)
00219 cflux[i] = cn[i]*bb[i];
00220 }
00221
00222
00223
00224
00225
00226
00227 static float fraction (float x, float y, float r_out) {
00228 float r,t,x_a,x_b,frac,tanao2,cosa,tanp2a,sqrt2o2;
00229
00230 r = sqrtf(x*x + y*y);
00231 sqrt2o2 = 0.5*M_SQRT2;
00232
00233
00234
00235 if(r > r_out+sqrt2o2)
00236 return(0.0);
00237
00238
00239
00240 if(r < r_out-sqrt2o2)
00241 return(1.0);
00242
00243
00244
00245
00246 x = fabsf(x);
00247 y = fabsf(y);
00248 if(y > x) {
00249 t = x;
00250 x = y;
00251 y = t;
00252 }
00253
00254
00255
00256 if (x > 0.0 && y > 0.0) {
00257 tanao2 = 0.5*y/x;
00258 tanp2a = x/y;
00259 cosa = x/sqrt(x*x + y*y);
00260 } else {
00261 tanao2 = 0.00005;
00262 tanp2a = 10000.0;
00263 cosa = 1.0;
00264 }
00265
00266
00267
00268 x_a = x - tanao2 + (r_out - r)/cosa;
00269 if(x_a < x+0.5) {
00270
00271
00272
00273 x_b = x + tanao2 + (r_out - r)/cosa;
00274
00275
00276
00277 if(x_a < x-0.5)
00278 frac = 0.5*MAX(0.0,x_b-(x-0.5))*MAX(0.0,x_b-(x-0.5))*tanp2a;
00279 else {
00280 if(x_b > x+0.5)
00281 frac = 1.0 - 0.5*(x+0.5-x_a)*(x+0.5-x_a)*tanp2a;
00282 else
00283 frac = 0.5-(x-x_a)+0.5*(x_b-x_a);
00284 }
00285 } else
00286 frac = 1.0;
00287
00288 return(frac);
00289 }
00290
00291
00292
00293 static void dchole (double a[IMNUM+1][IMNUM+1], double b[IMNUM+1], int n) {
00294 double sum, l[IMNUM+1][IMNUM+1], y[IMNUM+1];
00295 double aveigv, offset;
00296 int i, j, k;
00297
00298 restart:
00299 l[0][0] = sqrt(a[0][0]);
00300
00301 for(k = 1; k < n; k++) {
00302 for(j = 0; j <= k-1; j++) {
00303 sum = a[j][k];
00304 if(j != 0)
00305 for(i = 0; i <= j-1; i++)
00306 sum -= l[i][k]*l[i][j];
00307 l[j][k] = sum/l[j][j];
00308 }
00309 sum = a[k][k];
00310 for(i = 0; i <= k-1; i++)
00311 sum -= l[i][k]*l[i][k];
00312 if(sum <= 0.0) {
00313
00314 aveigv = a[0][0];
00315 for(i = 1; i < n; i++)
00316 aveigv += a[i][i];
00317
00318 offset = 0.1*aveigv/((double) n);
00319 for(i = 0; i < n; i++)
00320 a[i][i] += offset;
00321
00322 goto restart;
00323 }
00324 l[k][k] = sqrt(sum);
00325 }
00326
00327
00328
00329 y[0] = b[0]/l[0][0];
00330 for(i = 1; i < n; i++) {
00331 sum = b[i];
00332 for(k = 0; k <= i-1; k++)
00333 sum -= l[k][i]*y[k];
00334 y[i] = sum/l[i][i];
00335 }
00336
00337
00338
00339 b[n-1] = y[n-1]/l[n-1][n-1];
00340 for(i = n-2; i >= 0; i--) {
00341 sum = y[i];
00342 for(k = i+1; k < n; k++)
00343 sum -= l[i][k]*b[k];
00344 b[i] = sum/l[i][i];
00345 }
00346 }
00347
00348
00349
00350
00351
00352
00353
00354
00355