00001
00002
00003
00004
00005
00006
00007 #include <stdio.h>
00008 #include <math.h>
00009 #include "imcore.h"
00010 #include "util.h"
00011 #include "floatmath.h"
00012
00013
00014
00015 #define COL_NUMBER 1
00016 #define COL_FLUXISO 2
00017 #define COL_FLUXTOTAL 3
00018 #define COL_APFLUX1 4
00019 #define COL_X 5
00020 #define COL_Y 6
00021 #define COL_SIGMA 7
00022 #define COL_ELLIPT 8
00023 #define COL_PA 9
00024 #define COL_PEAKHEIGHT 10
00025 #define COL_AREAL1 11
00026 #define COL_AREAL2 12
00027 #define COL_AREAL3 13
00028 #define COL_AREAL4 14
00029 #define COL_AREAL5 15
00030 #define COL_AREAL6 16
00031 #define COL_AREAL7 17
00032 #define COL_AREAL8 18
00033 #define COL_APFLUX2 19
00034 #define COL_APFLUX3 20
00035 #define COL_APFLUX4 21
00036 #define COL_APFLUX5 22
00037 #define COL_RA 23
00038 #define COL_DEC 24
00039 #define COL_CLASS 25
00040 #define COL_STAT 26
00041 #define COL_APFLUX6 27
00042 #define COL_SKYLEV 28
00043 #define COL_SKYRMS 29
00044
00045
00046
00047 #define NCOLS 32
00048
00049
00050
00051 #define NRADS 6
00052 static float rmults[] = {0.5,1.0,M_SQRT2,2.0,2.0*M_SQRT2,4.0};
00053 static int nrcore = 1;
00054 static float apertures[NRADS];
00055
00056
00057
00058 static const char *ttype[NCOLS]={"No.","Isophotal_flux","Total_flux","Core_flux",
00059 "X_coordinate","Y_coordinate","Gaussian_sigma",
00060 "Ellipticity","Position_angle","Peak_height",
00061 "Areal_1_profile","Areal_2_profile","Areal_3_profile",
00062 "Areal_4_profile","Areal_5_profile","Areal_6_profile",
00063 "Areal_7_profile","Areal_8_profile","Core1_flux",
00064 "Core2_flux","Core3_flux","Core4_flux",
00065 "RA","DEC","Classification","Statistic",
00066 "Core5_flux","Skylev",
00067 "Skyrms","Blank_30","Blank_31","Blank_32"};
00068 static const char *tunit[NCOLS]={" ","Counts","Counts","Counts","Pixels","Pixels",
00069 "Pixels"," ","Degrees","Counts","Pixels","Pixels",
00070 "Pixels","Pixels","Pixels","Pixels","Pixels","Pixels",
00071 "Counts","Counts","Counts","Counts",
00072 "Radians","Radians","Flag","N-sigma",
00073 "Counts","Counts","Counts","Blank_30",
00074 "Blank_31","Blank_32"};
00075 static cpl_type tform[NCOLS]={CPL_TYPE_INT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00076 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00077 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00078 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00079 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00080 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00081 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00082 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00083 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00084 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,CPL_TYPE_FLOAT,
00085 CPL_TYPE_FLOAT,CPL_TYPE_FLOAT};
00086
00087 static int areal_cols[NAREAL] = {COL_AREAL1,COL_AREAL2,COL_AREAL3,COL_AREAL4,
00088 COL_AREAL5,COL_AREAL6,COL_AREAL7,COL_AREAL8};
00089
00090 extern void tabinit_1(void) {
00091
00092
00093
00094 tabinit_gen(NCOLS,ttype,tunit,tform);
00095
00096
00097
00098 imcore_xcol = COL_X;
00099 imcore_ycol = COL_Y;
00100 }
00101
00102 extern int do_seeing_1(ap_t *ap) {
00103 int retval,i;
00104 char *areal_colnames[NAREAL];
00105
00106
00107
00108 for (i = 0; i < NAREAL; i++)
00109 areal_colnames[i] = (char *)ttype[areal_cols[i]-1];
00110
00111
00112
00113 retval = do_seeing_gen(ap,ttype[COL_ELLIPT-1],ttype[COL_PEAKHEIGHT-1],
00114 areal_colnames);
00115
00116
00117
00118 return(retval);
00119 }
00120
00121 extern int process_results_1(ap_t *ap) {
00122 float momresults[8],ttotal,parmall[IMNUM][NPAR],ratio,cflux[NRADS*IMNUM];
00123 float sxx,syy,srr,sxy,ecc,temp,xx,theta,radeg,ell,iso_flux,total_flux;
00124 float apflux1,apflux2,apflux3,apflux4,apflux5,yy,sigma,peak,areal1,apflux6;
00125 float areal2,areal3,areal4,areal5,areal6,areal7,areal8,aa,bb;
00126 float skylev,skyrms,badpix[IMNUM];
00127 int iareal[NAREAL],nbit,i,k,nn,nr,mbit,j;
00128 long nrows;
00129 plstruct *pl;
00130 unsigned char *mflag;
00131
00132
00133
00134 moments(ap,momresults);
00135 if (momresults[0] < 0)
00136 return(VIR_FATAL);
00137 areals(ap,iareal);
00138
00139
00140
00141
00142 if (iareal[0] < ap->ipnop || momresults[3] < ap->xintmin)
00143 return(VIR_OK);
00144
00145
00146
00147 extend(ap,momresults[3],momresults[1],momresults[2],
00148 momresults[4],momresults[5],momresults[6],
00149 (float)iareal[0],momresults[7],&ttotal);
00150 ratio = MAX(ttotal,momresults[3])/momresults[3];
00151
00152
00153
00154 if (iareal[0] >= ap->mulpix && ap->icrowd)
00155 overlp(ap,parmall,&nbit,momresults[1],momresults[2],
00156 momresults[3],iareal[0],momresults[7]);
00157 else
00158 nbit = 1;
00159 if (nbit == 1) {
00160 parmall[0][0] = momresults[3];
00161 parmall[0][1] = momresults[1];
00162 parmall[0][2] = momresults[2];
00163 parmall[0][3] = ap->thresh;
00164 for (i = 4; i < 8; i++)
00165 parmall[0][i] = momresults[i];
00166 for (i = 0; i < NAREAL; i++)
00167 parmall[0][i+8] = (float)iareal[i];
00168 } else {
00169 mbit = 0;
00170 for (i = 0; i < nbit; i++) {
00171 if (parmall[i][1] > 1.0 && parmall[i][1] < ap->lsiz &&
00172 parmall[i][2] > 1.0 && parmall[i][2] < ap->csiz) {
00173 for (j = 0; j < NPAR; j++)
00174 parmall[mbit][j] = parmall[i][j];
00175 mbit++;
00176 }
00177 }
00178 nbit = mbit;
00179 }
00180
00181
00182
00183 for (i = 0; i < NRADS; i++)
00184 apertures[i] = rmults[i]*(ap->rcore);
00185
00186
00187
00188 for (i = 0; i < nbit; i++)
00189 badpix[i] = 0.0;
00190
00191
00192
00193 phopt(ap,parmall,nbit,NRADS,apertures,cflux,badpix,nrcore);
00194
00195
00196
00197 radeg = 180.0/M_PI;
00198 for (k = 0; k < nbit; k++) {
00199 sxx = parmall[k][4];
00200 sxy = parmall[k][5];
00201 syy = parmall[k][6];
00202 if(sxy > 0.0)
00203 sxy = MAX(1.0e-4,MIN(sxy,sqrtf(sxx*syy)));
00204 else
00205 sxy = MIN(-1.0e-4,MAX(sxy,-sqrtf(sxx*syy)));
00206
00207 srr = MAX(0.5,sxx+syy);
00208 ecc = sqrtf((syy-sxx)*(syy-sxx)+4.0*sxy*sxy)/srr;
00209 temp = MAX((1.0-ecc)/(1.0+ecc),0.0);
00210 ell = 1.0 - sqrtf(temp);
00211 ell = MIN(0.99,MAX(0.0,ell));
00212 xx = 0.5*(1.0+ecc)*srr-sxx;
00213 if(xx == 0.0)
00214 theta = 0.0;
00215 else
00216 theta = 90.0-radeg*atanf(sxy/xx);
00217
00218
00219
00220 nrows = cpl_table_get_nrow(tab);
00221 nobjects++;
00222 if (nobjects > nrows)
00223 (void)cpl_table_set_size(tab,nrows+INITROWS);
00224 nr = nobjects - 1;
00225 iso_flux = parmall[k][0];
00226 total_flux = ratio*parmall[k][0];
00227 apflux1 = cflux[k*NRADS + 0];
00228 apflux2 = cflux[k*NRADS + 1];
00229 apflux3 = cflux[k*NRADS + 2];
00230 apflux4 = cflux[k*NRADS + 3];
00231 apflux5 = cflux[k*NRADS + 4];
00232 apflux6 = cflux[k*NRADS + 5];
00233 xx = parmall[k][1];
00234 yy = parmall[k][2];
00235 sigma = sqrt(srr);
00236 peak = parmall[k][7];
00237 areal1 = parmall[k][8];
00238 areal2 = parmall[k][9];
00239 areal3 = parmall[k][10];
00240 areal4 = parmall[k][11];
00241 areal5 = parmall[k][12];
00242 areal6 = parmall[k][13];
00243 areal7 = parmall[k][14];
00244 if (nbit > 1 && k == 0)
00245 areal8 = 0.0;
00246 else
00247 areal8 = parmall[k][15];
00248 imcore_backest(ap,xx,yy,&skylev,&skyrms);
00249
00250
00251
00252 cpl_table_set_int(tab,ttype[COL_NUMBER-1],nr,nobjects);
00253 cpl_table_set_float(tab,ttype[COL_FLUXISO-1],nr,iso_flux);
00254 cpl_table_set_float(tab,ttype[COL_FLUXTOTAL-1],nr,total_flux);
00255 cpl_table_set_float(tab,ttype[COL_APFLUX1-1],nr,apflux2);
00256 cpl_table_set_float(tab,ttype[COL_X-1],nr,xx);
00257 cpl_table_set_float(tab,ttype[COL_Y-1],nr,yy);
00258 cpl_table_set_float(tab,ttype[COL_SIGMA-1],nr,sigma);
00259 cpl_table_set_float(tab,ttype[COL_ELLIPT-1],nr,ell);
00260 cpl_table_set_float(tab,ttype[COL_PA-1],nr,theta);
00261 cpl_table_set_float(tab,ttype[COL_PEAKHEIGHT-1],nr,peak);
00262 cpl_table_set_float(tab,ttype[COL_AREAL1-1],nr,areal1);
00263 cpl_table_set_float(tab,ttype[COL_AREAL2-1],nr,areal2);
00264 cpl_table_set_float(tab,ttype[COL_AREAL3-1],nr,areal3);
00265 cpl_table_set_float(tab,ttype[COL_AREAL4-1],nr,areal4);
00266 cpl_table_set_float(tab,ttype[COL_AREAL5-1],nr,areal5);
00267 cpl_table_set_float(tab,ttype[COL_AREAL6-1],nr,areal6);
00268 cpl_table_set_float(tab,ttype[COL_AREAL7-1],nr,areal7);
00269 cpl_table_set_float(tab,ttype[COL_AREAL8-1],nr,areal8);
00270 cpl_table_set_float(tab,ttype[COL_APFLUX2-1],nr,apflux1);
00271 cpl_table_set_float(tab,ttype[COL_APFLUX3-1],nr,apflux3);
00272 cpl_table_set_float(tab,ttype[COL_APFLUX4-1],nr,apflux4);
00273 cpl_table_set_float(tab,ttype[COL_APFLUX5-1],nr,apflux5);
00274 cpl_table_set_float(tab,ttype[COL_APFLUX6-1],nr,apflux6);
00275 cpl_table_set_float(tab,ttype[COL_SKYLEV-1],nr,skylev);
00276 cpl_table_set_float(tab,ttype[COL_SKYRMS-1],nr,skyrms);
00277
00278
00279
00280 if (ellfp != NULL) {
00281 aa = sqrtf(areal1/(M_PI*(1.0 - ell)));
00282 bb = aa*(1.0 - ell);
00283 fprintf(ellfp,"image; ellipse %9.3f %9.3f %7.2f %7.2f %6.1f\n",
00284 xx,yy,aa,bb,theta);
00285 }
00286 }
00287
00288
00289
00290 pl = ap->plarray;
00291 mflag = ap->mflag;
00292 for (k = 0; k < ap->npl_pix; k++) {
00293 nn = ap->lsiz*(pl[k].y - 1) + pl[k].x - 1;
00294 mflag[nn] = MF_POBJPIX;
00295 }
00296
00297
00298
00299 return(VIR_OK);
00300 }
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338