#include <stdio.h>
#include <math.h>
#include "export.h"

/* To compile, link, and use:
     cc -K pic -G -I /opt/local/rsi/idl/external -c cgetrng.c
     cc -G -o cgetrng cgetrng.o

     linkimage,'CGETRNG','/gryll/data1/buie/idl/Custom/cgetrng',0,'cgetrng', $
        max_args=8,min_args=8

     After this statement is executed, cgetrng will be drawn from the custom
     routine and not the IDL cgetrng.pro file.
*/

float cg_rnd(x) /* round a floating point number to the nearest whole number */

float x;

{
   if ( x < 0 ) return( (int)(x-0.5) );
   else return( (int)(x+0.5) );
}

/* This procedure is called to determine how to iterate when integrating
 * over a circle.  The circle's center is at (xc,yc), and its radius is r.
 * For pixels with x-coordinate x, those in the intervals [y0,y1) and [y2,y3)
 * are on or near the circle.  Those in the interval [y1,y2) are definitely
 * inside; all others are definitely outside.
 *    Of course, the routine can be called to determine an interval for fixed
 * y by calling it as cgetrng(yc,xc,r,y,x0,x1,x2,x3).
 *    The appropriate way to integrate over a circle is therefore as follows:
 * cgetrng(xc,yc,r,round(xc),y0,y1,y2,y3);
 * for (y = y0; y <= y3-1; y=y+1) {
 *    cgetrng (yc, xc, r, y, x0, x1, x2, x3);
 *    for (x = x0; x <= x1-1; x=x+1) sum = sum + value(x,y)*pixwt(xc,yc,r,x,y);
 *    for (x = x1; x <= x2-1; x=x+1) sum = sum + value(x,y);
 *    for (x = x2; x <= x3-1; x=x+1) sum = sum + value(x,y)*pixwt(xc,yc,r,x,y);
 *    }
 */

#define sqrt2   1.414213563     /* rounded up to increase size of
                                   uncertain areas  */

void cgetrng(int argc, IDL_VPTR argv[])

/* cgetrng(xc, yc, r, x, y0, y1, y2, y3) */

/* float   xc, yc, r;              Center and radius of the circle. */
/* int     x;                      X coordinate of the pixels. */
/* int     *y0, *y1, *y2, *y3;     Boundaries of circle in Y. */

{
   float   a,b,outdsq,indsq,outd,ind;

   IDL_VPTR vxc,vyc,vr,vx,vy0,vy1,vy2,vy3;
   float xc,yc,r,x;
   long  y0,y1,y2,y3;

   vxc = IDL_CvtFlt(1,&argv[0]);
   vyc = IDL_CvtFlt(1,&argv[1]);
   vr  = IDL_CvtFlt(1,&argv[2]);
   vx  = IDL_CvtFlt(1,&argv[3]);

   xc = vxc->value.f;
   yc = vyc->value.f;
   r  = vr->value.f;
   x  = vx->value.f;

   if (r <= 0) /* then it misses completely */ outdsq = -1;
   else {
      a = r*r + 0.5 - (x-xc)*(x-xc);
      b = sqrt2 * r;
      outdsq = a + b;
      if(b < 1) /*then indsq would be invalid--say no interior*/ indsq= -1;
      else /* formula works */ indsq = a-b; }

   if (outdsq < 0) /* complete miss */ y0 = y1 = y2 = y3 = cg_rnd(yc);
   else /* there is some intersection */ {
      outd = sqrt(outdsq);
      y0 = ceil(yc-outd);
      y3 = floor(yc+outd) + 1;
      if (indsq < 0) /* no interior */ y1 = y2 = cg_rnd(yc);
      else /* there is a certain interior */ {
         ind = sqrt(indsq);
         y1 = ceil(yc-ind);
         y2 = floor(yc+ind) + 1;}
      }

   vy0 = IDL_Gettmp();
   vy1 = IDL_Gettmp();
   vy2 = IDL_Gettmp();
   vy3 = IDL_Gettmp();

   vy0->type=IDL_TYP_LONG;
   vy1->type=IDL_TYP_LONG;
   vy2->type=IDL_TYP_LONG;
   vy3->type=IDL_TYP_LONG;

   vy0->value.l = y0;
   vy1->value.l = y1;
   vy2->value.l = y2;
   vy3->value.l = y3;

   IDL_VarCopy(vy0,argv[4]);
   IDL_VarCopy(vy1,argv[5]);
   IDL_VarCopy(vy2,argv[6]);
   IDL_VarCopy(vy3,argv[7]);

   IDL_DELTMP(vxc);
   IDL_DELTMP(vyc);
   IDL_DELTMP(vr);
   IDL_DELTMP(vx);
}