← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to prim/general/libsrc/cjmagn.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
prim/general/libsrc/cjmagn.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1995-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
        Internet e-mail: midas@eso.org
 
21
        Postal address: European Southern Observatory
 
22
                        Data Management Division 
 
23
                        Karl-Schwarzschild-Strasse 2
 
24
                        D 85748 Garching bei Muenchen 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
.IDENTIFIER  Cjmagn
 
30
.LANGUAGE    C
 
31
.AUTHOR      K. Banse           ESO - Garching
 
32
.KEYWORDS    
 
33
.PURPOSE     calculate the magnitude of a star using method JMETH
 
34
.ALGORITHM   
 
35
             
 
36
.IN/OUTPUT
 
37
  call as    stat = Cjmagn( jmeth, arr, npix, ni, nb, xycen,
 
38
                            mag, dmag, sky, dsky, apix, flux )
 
39
  input:
 
40
          int   jmeth           method of magnitude calculation
 
41
                                 1: sum of 9 central pixels
 
42
                                 2: sum of all pixels in "flux" area (DEFAULT)
 
43
                                 3: circular aperture
 
44
          float *arr            data around the star
 
45
          int   npix[2]         dimension of ARR
 
46
          int   ni              number of pixels between "flux" area 
 
47
                                and background
 
48
          int   nb              width of background along edge
 
49
          float xycen[2]        center of star in pixels in array "arr" (x,y)
 
50
  output:
 
51
          float *mag            magnitude of star
 
52
          float *dmag           uncertainty in magnitude determination
 
53
          float *sky            mean sky intensity
 
54
          float *dsky           uncertainty in sky determination
 
55
          float *nrpix          number of pixels used for magn. determination
 
56
          float *flux           total flux of the star
 
57
 
 
58
.RETURNS     status      1: put star in center of image
 
59
                         0: ok
 
60
                        -1: array "arr" too small
 
61
                        -2: center not in array "arr"
 
62
.ENVIRONment MIDAS
 
63
             #include <midas_def.h>      Prototypes for MIDAS interfaces
 
64
 
 
65
.VERSIONS    1.00       940411  from JMAGN.FOR    RvH
 
66
 
 
67
 090703         last modif
 
68
------------------------------------------------------------*/
 
69
 
 
70
/* Define _POSIX_SOURCE to indicate that this is a POSIX program */
 
71
 
 
72
#define  _POSIX_SOURCE 1
 
73
 
 
74
 
 
75
#include <math.h>
 
76
 
 
77
/*
 
78
 
 
79
*/
 
80
 
 
81
/*++++++++++++++++++++++++++++++
 
82
.IDENTIFIER  static function pixfrac
 
83
.PURPOSE     calculate the fraction of a pixel which is inside a circle
 
84
             with radius r around (xc,yc)
 
85
.IN/OUTPUT
 
86
  call as    static function
 
87
  input: 
 
88
          double f              pixel value (flux)
 
89
          double df[4]          neg. partial derivative of f with respect to x
 
90
                                pos. partial derivative of f with respect to x
 
91
                                neg. partial derivative of f with respect to y
 
92
                                pos. partial derivative of f with respect to y
 
93
          float  xdelt          xdelta (x - ycenter)
 
94
          float  ydelt          ydelta (y - ycenter)
 
95
          double r              radius of circle
 
96
  output:
 
97
          double *fi            flux inside circle
 
98
          double *pi            area inside circle
 
99
------------------------------*/
 
100
 
 
101
#ifdef __STDC__
 
102
      static void PFRAC(double f, double df[4], float xdelt, float ydelt,
 
103
                        double r, double *fi, double *pi)
 
104
#else
 
105
      static void PFRAC(f, df, xdelt, ydelt, r, fi, pi)
 
106
      float  xdelt, ydelt;
 
107
      double f, df[4], r, *fi, *pi;
 
108
#endif
 
109
 
 
110
{
 
111
double fx, fy, x, y, ai;
 
112
double v, vy, rx, ry, rp;
 
113
double dx = 0.1, dy = 0.1, da = 0.01;
 
114
 
 
115
 
 
116
*fi = *pi = 0.0;
 
117
ai = f - 0.5 * (df[0] - df[1] + df[2] - df[3]);
 
118
y  = -0.5 + (0.5 * dy);
 
119
 
 
120
 
 
121
/* Subdivide pixel and check */
 
122
 
 
123
sect_100:
 
124
if (y > 0.0)
 
125
   fy = df[2];
 
126
else
 
127
   fy = df[3];
 
128
x = -0.5 + (0.5 * dx);
 
129
vy = fy * y;
 
130
ry = (ydelt + y);
 
131
ry *= ry;
 
132
 
 
133
sect_200:
 
134
if (x > 0.0) 
 
135
   fx = df[0];
 
136
else
 
137
   fx = df[1];
 
138
v = vy + (fx * x);
 
139
rx = (xdelt + x);
 
140
rx *= rx; 
 
141
rp = sqrt(rx + ry);
 
142
 
 
143
 
 
144
/* check if inside */
 
145
 
 
146
if ((r - rp) >= 0.0)
 
147
   {
 
148
   *fi += v;
 
149
   *pi += da;
 
150
   }
 
151
 
 
152
x += dx;
 
153
if (x < 0.5) goto sect_200;
 
154
 
 
155
y += dy;
 
156
if (y < 0.5) goto sect_100;
 
157
 
 
158
 
 
159
*fi = (ai * *pi) + (*fi * da);
 
160
}
 
161
 
 
162
/*
 
163
 
 
164
*/
 
165
 
 
166
int Cjmagn( jmeth, arr, npix, ni, nb, pfac, xycen,
 
167
            mag, dmag, sky, dsky, nrpix, flux )
 
168
int   jmeth, *npix, ni, nb;
 
169
float *arr, *pfac, *xycen, *mag, *dmag, *sky, *dsky, *nrpix, *flux;
 
170
 
 
171
{
 
172
register int ii, ix, iy;
 
173
int   nn, ninb, nxb[2], nyb[2];
 
174
int   ierr, xmin, ymin, nrx, nry;
 
175
 
 
176
float  rn, rns, dflux;
 
177
float  *pntra, *pntr, a, xrval, yrval, fac;
 
178
 
 
179
double s, sk, sq, akap;
 
180
 
 
181
 
 
182
rns = 0.0;
 
183
fac = *pfac;
 
184
*mag  = -9999;                                  /* set return values */
 
185
*dmag = *sky = *dsky = *flux = 0.0;
 
186
*nrpix = 0.0;
 
187
ierr = 0;
 
188
ninb = ni + nb;
 
189
nn = (2 * ninb) + 3;
 
190
 
 
191
 
 
192
/* check dimensions and position */
 
193
 
 
194
if ((npix[0] <nn) || (npix[1] < nn)) return (-3);        /* array too small */
 
195
 
 
196
if ((*xycen < (ninb-1)) || (xycen[0] > (float)(npix[0]-ninb-1)))
 
197
   {
 
198
   if (jmeth == 1) return (-2);                 /* center not in array */
 
199
   ierr = 1;                                    /* put star in center */
 
200
   xycen[0] = (npix[0]-1) * 0.5;
 
201
   }
 
202
 
 
203
if ((xycen[1] < (ninb-1)) || (xycen[1] > (float)(npix[1]-ninb-1))) 
 
204
   {
 
205
   if (jmeth == 1) return (-2);                 /* center not in array */
 
206
   ierr = 1;                                    /* put star in center */
 
207
   xycen[1] = (npix[1]-1) * 0.5;
 
208
   }
 
209
  
 
210
 
 
211
/* for all methods we determine the sky intensity by clipping
 
212
   with 2.0 * sigma */
 
213
 
 
214
akap = 1.0e+30;                         /* huge Kappa */
 
215
nxb[0] = nb - 1;
 
216
nxb[1] = npix[0] - nb;
 
217
nyb[0] = nb - 1;
 
218
nyb[1] = npix[1] - nb;
 
219
 
 
220
 
 
221
/* calculate background sky, if so required */
 
222
 
 
223
if (nb > 0)     
 
224
   {
 
225
   xrval = 0.0;                         /* xrval = 0.0 (initial *sky) */
 
226
 
 
227
   for (ii=0; ii<10; ii++)                      /* 10 iterations... */
 
228
      {
 
229
      s = 0.0;
 
230
      sq = 0.0;
 
231
      nn = 0;
 
232
      pntr = arr;
 
233
 
 
234
      for (ix=0; ix<npix[0]; ix++)
 
235
         {
 
236
         for (iy=0; iy<npix[1]; iy++)
 
237
            {
 
238
            if ((ix <= nxb[0]) || (ix >= nxb[1]) ||
 
239
                (iy <= nyb[0]) || (iy >= nyb[1]))
 
240
               {
 
241
               a = *pntr;
 
242
               if ( fabs((double) (a - xrval) ) <= akap )
 
243
                  {
 
244
                  nn++;
 
245
                  s  += a;
 
246
                  sq += (a*a);
 
247
                  }
 
248
               }
 
249
            pntr++;
 
250
            }
 
251
         }
 
252
 
 
253
      if (nn > 0)
 
254
         {
 
255
         rns = (float) (nn);
 
256
         xrval  = (float) (s / rns);
 
257
         a = sq/rns - (xrval*xrval);
 
258
         if (a > 0.0)
 
259
            {
 
260
            *dsky = (float) sqrt((double) a); 
 
261
            akap  = fac * (*dsky);
 
262
            }
 
263
         else
 
264
            {
 
265
            akap = *dsky = 0.0; 
 
266
            }
 
267
         }
 
268
      }
 
269
 
 
270
   *sky  = xrval;
 
271
   }
 
272
 
 
273
 
 
274
/* calculate the magnitude for a given method */
 
275
 
 
276
if (jmeth != 3) 
 
277
   {
 
278
 
 
279
   /* method #1: sum of 9 central pixels */
 
280
   
 
281
   if (jmeth == 1)
 
282
      {
 
283
      xmin = (int) floor((double) xycen[0] );   /* C indexing, so remove */
 
284
      nrx = xmin + 2;                   /* subtraction of 1.0 from xycen */
 
285
      ymin = (int) floor((double) xycen[1] );
 
286
      nry = ymin + 2;
 
287
      }
 
288
 
 
289
   /* method #2: sum of all pixels within background */
 
290
 
 
291
   else   
 
292
      {
 
293
      xmin = ymin = ninb + 1;
 
294
      nrx = npix[0] - ninb;
 
295
      nry = npix[1] - ninb;
 
296
      }
 
297
 
 
298
   /* sum flux within area */
 
299
 
 
300
   nn = 0;
 
301
   s = 0.0;
 
302
 
 
303
   if ((xmin > nrx) || (ymin > nry)) return (-3);
 
304
 
 
305
   pntra = arr + (xmin-1 + (npix[0]*(ymin-1)));         /* point to start */
 
306
   for (iy=0; iy<(nry-ymin+1); iy++)
 
307
      {
 
308
      pntr = pntra;
 
309
      for (ix=0; ix<(nrx-xmin+1); ix++) 
 
310
         {
 
311
         nn ++;
 
312
         s += *pntr++;
 
313
         }
 
314
      pntra += npix[0];
 
315
      }
 
316
 
 
317
   rn = (float)(nn);
 
318
   }
 
319
 
 
320
 
 
321
/* method #3: circular aperture */
 
322
 
 
323
else
 
324
   {
 
325
   int      npxm1, npym1;
 
326
   float    ry;
 
327
   double   fi, r, dr, ro, ri, pi, val, dval[4];
 
328
 
 
329
   rns = 0.0;
 
330
   ro = npix[0] - xycen[0] - 1.0 - nb;
 
331
   ri = ro - ni;
 
332
   dr = sqrt((double)0.5);
 
333
   sk = 0.0;
 
334
   sq = 0.0;
 
335
   s = 0.0;
 
336
   rn = 0.0;
 
337
   npxm1 = npix[0] - 1;
 
338
   npym1 = npix[1] - 1;
 
339
 
 
340
   pntra = arr;
 
341
   for (iy=0; iy<npix[1]; iy++)
 
342
      {
 
343
      pntr = pntra;
 
344
      yrval = (float) (iy - xycen[1]);
 
345
      ry = yrval * yrval;
 
346
      for (ix=0; ix<npix[0]; ix++)
 
347
         {
 
348
         val = *pntr;
 
349
         xrval = (float) (ix-xycen[0]);
 
350
         r = sqrt((double) ((xrval*xrval) + ry));
 
351
         if ( (nb > 0) && (r >= ro) && (fabs(val - *sky) <= akap) )     
 
352
            {
 
353
            rns += 1.0;
 
354
            sk += val;
 
355
            sq += (val*val);
 
356
            }
 
357
               
 
358
         if ((r - dr) < ri) 
 
359
            {
 
360
            if ((r + dr) <= ri) 
 
361
               {
 
362
               s += val;
 
363
               rn += 1.0;
 
364
               }
 
365
            else
 
366
               {
 
367
               dval[0] = *(pntr+1) - val;
 
368
               a = val - *(pntr-1);
 
369
               if (ix == npxm1) dval[0] = a;
 
370
               dval[1] = a;
 
371
               if (ix == 0) dval[1] = dval[0];
 
372
 
 
373
               dval[2] = *(pntr+npix[0]) - val;
 
374
               a = val - *(pntr-npix[0]);
 
375
               if (iy == npym1) dval[2] = a;
 
376
               dval[3] = a;
 
377
               if (iy == 0) dval[3] = dval[2];
 
378
 
 
379
               PFRAC(val,dval,xrval, yrval,ri,&fi,&pi);
 
380
               s  += fi;
 
381
               rn += pi;
 
382
               }
 
383
            }
 
384
         pntr++;
 
385
         }
 
386
             
 
387
      pntra += npix[0];
 
388
      }
 
389
      
 
390
   if (rns > 1.0)
 
391
      { 
 
392
      xrval = (float) (sk / rns);
 
393
      *sky = xrval;
 
394
      a = (float) (sq/rns - (xrval*xrval));
 
395
      if (a > 0.0)
 
396
         *dsky = (float) sqrt((double) a);
 
397
      else
 
398
         *dsky = 0.0; 
 
399
      }
 
400
   }
 
401
      
 
402
 
 
403
/* For all methods: compute magnitude and uncertainty */
 
404
 
 
405
*nrpix = rn;
 
406
 
 
407
if (nb > 0)
 
408
   {
 
409
   *flux = s - (rn * *sky);
 
410
   dflux = rn * (*dsky) * (float) sqrt((double) (1.0/rn + 1.0/rns) );
 
411
   }
 
412
else
 
413
   {
 
414
   *flux = s;
 
415
   dflux = 0.;
 
416
   }
 
417
 
 
418
if (*flux >= (0.1*dflux))     
 
419
   {
 
420
   *mag  = -2.5 * log10( (double) *flux );
 
421
   *dmag = 1.0857362 * dflux / *flux;
 
422
   }
 
423
else
 
424
   ierr = 2;
 
425
 
 
426
return ierr;
 
427
}
 
428