//////////////////////////////////////////////////////////////////////////
//
//	test stuff for TMO_Transit
//
/////////////////////////////////////////////////////////////////////////


#include "SpatialIndex.h"
#include "SpatialEdge.h"
#include "SpatialGeneral.h"
#include "SpatialVector.h"
#include "SpatialConstraint.H"
#include "SpatialConvex.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NORTH 1
#define SOUTH 0
#define INDEX_SIZE 14

static int adata4[600][1280], bdata4[600][1280], triangle[10000];
static int cdata4[600][1280], ddata4[600][1280], edata4[600][1280];
static short adata2[600][1280], bdata2[600][1280], cdata2[600][1280];
static short ddata2[600][1280], edata2[600][1280];
static float image[1600][400], imagecnt[1600][400];

typedef struct {
  float64 *response_;
  size_t length_;
  size_t capacity_;
} ResponseList;

typedef struct {
  int *hits_;
  size_t length_;
  size_t capacity_;
} HitList;

//////////////////////////////////////////////////////////////////////
//  Initialization for the lists
//////////////////////////////////////////////////////////////////////
void responseListInit( ResponseList *r ) 
{
  r->response_ = NULL;
  r->length_ = 0;
  r->capacity_ = 0;
}
void hitListInit( HitList *h ) 
{
  h->hits_ = NULL;
  h->length_ = 0;
  h->capacity_ = 0;
}
//////////////////////////////////////////////////////////////////////
//  appendation for the lists
//////////////////////////////////////////////////////////////////////
void responseListAppend( ResponseList *r, float64 i ) 
{
  size_t newsize = r->capacity_ * 2;
  
  if(r->length_ == r->capacity_) 
  {
    if(newsize == 0)newsize = 2;
    if(r->response_ == NULL)
      r->response_ = (float64 *)malloc( newsize * sizeof(float64) );
    else {
      r->response_ = (float64 *)realloc( r->response_, newsize * sizeof(float64) );
    }
    r->capacity_ = newsize;
  }
  r->response_[r->length_] = i;
  r->length_++;

  return;
}
void hitListAppend( HitList *h, int i ) 
{
  size_t newsize = h->capacity_ * 2;

  if(h->length_ == h->capacity_) 
  {
    if(newsize == 0)newsize = 2;
    if(h->hits_ == NULL)
      h->hits_ = (int *)malloc( newsize * sizeof(int) );
    else {
      h->hits_ = (int *)realloc( h->hits_, newsize * sizeof(int) );
    }
    h->capacity_ = newsize;
  }
  h->hits_[h->length_] = i;
  h->length_++;

  return;
}
//////////////////////////////////////////////////////////////////////
//  memory clearing for the lists
//////////////////////////////////////////////////////////////////////
void responseListClear( ResponseList *r , bool keep ) 
{
  if(!keep)
  {
    free(r->response_);
    r->capacity_ = 0;
  }
  r->length_ = 0;
}
void hitListClear( HitList *h , bool keep ) 
{
  if(!keep)
  {
    free(h->hits_);
    h->capacity_ = 0;
  }
  h->length_ = 0;
}
//////////////////////////////////////////////////////////////////////
//  calculate the center of a given spatial triangle
//////////////////////////////////////////////////////////////////////
void WspatialTriangleCenter( SpatialIndex *si, 
				IDTYPE *nodename, 
				float64 *ra, float64 *dec )
{
  SpatialVector *v1, *v2, *v3;
  IDTYPE bitname;
  v1 = spatialVectorNew();
  v2 = spatialVectorNew();
  v3 = spatialVectorNew();
  bitname = spatialIndexLeafNumberById(si, *nodename);
 
  //get vertices of triangle
  spatialIndexNodeVertexLeaf( si, bitname, v1, v2, v3 );

  //average over ra/dec of vertices to find center ra/dec
  *ra = ( v1->ra_ + v2->ra_ + v3->ra_ ) / 3;
  *dec = ( v1->dec_ + v2->dec_ + v3->dec_ ) / 3;
  
  if ( (abs(*ra - v1->ra_)) > 10 )
  {
    if (v1->ra_>180)  v1->ra_ = v1->ra_ - 360;      
    if (v2->ra_>180)  v2->ra_ = v2->ra_ - 360;      
    if (v3->ra_>180)  v3->ra_ = v3->ra_ - 360;
  
  *ra = ( v1->ra_ + v2->ra_ + v3->ra_ ) / 3;
  // 0<ra<360
  if(*ra<0)  *ra = *ra + 360;
  }
  
  spatialVectorDelete(v1);
  spatialVectorDelete(v2);
  spatialVectorDelete(v3);
  
  return;
}
//////////////////////////////////////////////////////////////////////
//	Euler rotation function
//////////////////////////////////////////////////////////////////////
void rotate(float64 alfa, float64 beta, float64 gamma, float64 *ra, float64 *dec)
{
  float64 sinb, cosb, sinl, cosl, sinp, cosp, cost;
  float64 degtorad=0.0174533, radtodeg=57.29578;
  float64 phi, rlat;
  
  phi = *ra;
  rlat = *dec;
  
  //convert from degrees to radians
  alfa = alfa * degtorad;
  beta = beta * degtorad;
  gamma = gamma * degtorad;
  phi = phi * degtorad;
  rlat = rlat * degtorad;
  
  phi = phi - alfa;  // rotation 1
  sinb = sin(beta);  // rotation 2...
  cosb = cos(beta);
  sinl = sin(rlat);
  cosl = cos(rlat);
  sinp = sin(phi);
  cosp = cos(phi);
  cost = ( sinl*cosb ) + ( cosl*sinb*cosp );
  
  if (abs(cost)!=1)
  {
    sinp = cosl*sinp;
    cosp = ( cosl*cosb*cosp ) - ( sinl*sinb );
  }else if (sinb = 0.0)
  {
    cosp = cosb*cosp;
  }
  
  phi = atan2(sinp,cosp);
  phi = phi - gamma;  // rotation 3
  
  //convert back to degrees
  alfa = alfa * radtodeg;
  beta = beta * radtodeg;
  gamma = gamma * radtodeg;
  phi = phi * radtodeg;
  rlat = rlat * radtodeg;
  
  phi = fmod(fmod(phi,360.0) + 360.0, 360.0);  // 0<phi<360
  rlat = 90.0 - ( acos(cost) * radtodeg );

  *ra = phi;
  *dec = rlat;

  return;
}
//////////////////////////////////////////////////////////////////////
// TEST: read in arbitrary file name, pedestal, print out file
//////////////////////////////////////////////////////////////////////

int mainwe(int argc,char ** argv)
{
  
  char nameFile[] = "TMO420000.nic";
  char nameoutFile[] = "Test_out_0000.grd";
  FILE  *inFile, *outFile;
  short head[3];
  size_t itemsread, newlen;
  float64 ped=0;
  int xarg, yarg, one=1;
  int pedcnt=0, i, j, k, l, h, ps=768000, *htemp;
  float64 aq, bq, cq, rq, thetax, thetaxsq, thetay, cx, cy, cz, crsq, cr, phi, rlat;
  float64 gamma=0.0, beta, alfa=0.50, camlat=34.9471, cra, cdec;
  float64 degperminute=0.2506844731,degtorad=0.0174533, radtodeg=57.29578;
  float64 x0=630.6, y0=1291.38, r0=1194.5, rccd, dr, ra0;
  float64 response=0.0, x[4], y[4], noderesponse, *pra, *pdec;
  float64 *rtemp;
  SpatialIndex *si;
  SpatialVector *v[4], *tempv;
  SpatialConvex *pixelconvex;
  SpatialMarkup *markup;
  SpatialIDVec *partialvec, *fullvec, *imagevec;
  HitList *hlist;
  ResponseList *rlist;
  IDTYPE *partialnode, *fullnode, *imagenode;
  
  
  printf("Starting...\n");

 for(k=170; k<171; k++)
 {
    sprintf(&nameFile[5],"%04d.nic",k);
    printf("inFile Name: %s\n",nameFile);
    
    sprintf(&nameoutFile[9],"%04d.grd",k);
    printf("outFile Name: %s\n",nameoutFile);
    
  // Read in both file
  inFile = fopen(&nameFile,"rb");
  itemsread = fread(&head, sizeof(short),3,inFile);  
  itemsread = fread(&adata2, sizeof(short),ps,inFile);
  fclose(inFile);


  for(i=0; i<600; i++)
  {
    for(j=0; j<1280; j++)
    {
      adata4[i][j] = 0;
    }
  }
  
  /*for(i=0; i<600; i++)
  {
    for(j=0; j<1280; j++)
    {
      adata4[i][j] = adata2[i][j];
      if(adata4[i][j]<0)
      {
        adata4[i][j] = adata4[i][j] + 65536;
      }
    }
  }*/
  
  
  for (i=0; i<600; i++)
  {
    for (j=0; j<50; j++)
    {                            
     rccd = sqrt((j-x0)*(j-x0) + (i-y0)*(i-y0));
     dr = fabs(rccd - r0);      
     //less than 1.5 deg x 18.0226 pix/deg
     if( dr < 27.0339 )
     {
      if( fabs((j-x0)/(i-y0)) < 0.57735 )// less than tan 30=0.57735
      {
	adata4[i][j] = adata2[i][j];
        if(adata4[i][j]<0)
        {
          adata4[i][j] = adata4[i][j] + 65536;
        }  
      }
     }
    }
   }
 
  /*pedestal subtraction from both files
  for(i=2; i<600; i++)
  {
    for(j=1271; j<1280; j++)
    {
      ped = ped + adata4[i][j];
      pedcnt++;
    }
    for(j=4; j<13; j++)
    {
      ped = ped + adata4[i][j];
      pedcnt++;
    }
  }  
  
  ped = ped/pedcnt;   
  printf("Ped: %f\n",ped);
  
  for(i=0; i<600; i++)
  {
    for(j=0; j<1280; j++)
    {
      adata4[i][j] = adata4[i][j] - ped;
    }
  }
  ped=0;
  pedcnt=0;
  */

  printf("Writing to file: 'Test_out.grd'\n");
  
  //open file
  outFile = fopen(nameoutFile,"w");
      
  //write header
  fprintf(outFile,"DSAA\n1280 600\n0 1280\n0 600\n0 65000\n");

  //write array
  for (i=0; i<600;)
  {
      for (j=0; j<1280;)
      {
         fprintf(outFile,"%8.6d ",adata4[i][j]);              
       j++;
      }
    i++;
  }
  fclose(outFile);
  
  
 }//end k-loop for each file name

  
  return 0;
}
//////////////////////////////////////////////////////////////////////
// TEST: work with coordinates, add precession
// Vega=[279.3,38.8], Enif=[326.1,9.875]
//////////////////////////////////////////////////////////////////////
int mainfg(int argc,char ** argv)
{
  int i, j, k, l, h, ps=768000, *htemp;
  float64 aq, bq, cq, rq, thetax, thetaxsq, thetay, cx, cy, cz, crsq, cr, phi, rlat;
  float64 gamma=0.0, beta, alfa=0.50, camlat=34.9471, cra, cdec;
  float64 degperminute=0.2506844731,degtorad=0.0174533, radtodeg=57.29578;
  float64 x0=630.6, y0=1291.38, r0=1194.5, rccd, ra0;
  float64 response=0.0, x[4], y[4], noderesponse, *pra, *pdec;
  float64 ratotal=0, dectotal=0, rafinal=0, decfinal=0;
  float64 am=1.2812323, an=0.556753, t;

          
       // adjustment for the night after 41xxx(42xxx)
       ra0=279.425;//**ad hoc**
       h=443;//frame number + 1
       
       //preession parameters from AQ4, p670
       t=0.017215;//time from J2000 in hundreds of y, for midnight 9/20/01
       am = am*t;
       an = an*t;
  
       //pixel coordinates of centroid
       i = 113;
       j = 783;
                                            
       x[0] = j - 0.5;
       y[0] = i - 0.5;
       x[1] = j - 0.5;
       y[1] = i + 0.5;
       x[2] = j + 0.5;
       y[2] = i + 0.5;
       x[3] = j + 0.5;
       y[3] = i - 0.5;
       
       //loop for each vertex of a pixel
       for (k=0; k<4; k++)
       {         
         rccd = sqrt((x[k] - x0)*(x[k] - x0) + (y[k] - y0)*(y[k] - y0));
         thetax = radtodeg*asin( (x[k] - x0) / rccd );
         thetax = thetax / 0.999744;//atmosphere adjustment
                        
         //quadratic coefficiants for the formula: rq=aqTY^2+bqTY+cq
         //from July 11, 2001 memo(ignoring 3rd and 4th order terms)
         /*
         thetaxsq = thetax*thetax;         
         aq = 0.1164 - ( pow(10,-5)*thetaxsq );
         bq = ( 3.6*pow(10,-5)*thetaxsq ) - 18.0266;
         cq = ( 0.000936*thetaxsq ) - ( 1.044*pow(10,-7)*thetaxsq*thetaxsq );
         */
         rq = (rccd - r0)/1.013;//y-plate scale adjustment  
         //thetay = ( -bq - sqrt( (bq*bq)-(4*aq*(cq-rq))) ) / ( 2*aq );
         thetay = rq / 18.0226; //note: +18.0226 to flip thetay/cy sign relative to Andy's program
                                                            
         cx = -sin(thetax*degtorad);
         cy = sin(thetay*degtorad);         
         crsq = (cx*cx)+(cy*cy);
         if(crsq>0)
         {
           phi = atan2(cy,cx);
         }else
         {
           phi = 0;
         }
         cz = sqrt(1-(crsq));
         rlat = radtodeg*asin(cz);
         phi = radtodeg*phi; 
                 
         beta = camlat-90;
         gamma = (-h*degperminute)-ra0; //(frame#+1)*degperminute  
         pra = &phi;
         pdec = &rlat;          
         
         rotate(alfa, beta, gamma, pra, pdec);  
         
         //printf("Pra, Pdec: %f %f\n",*pra,*pdec);
                 
         //precession adjustment
         *pra = *pra + 0.5*( am+an*sin(degtorad*(*pra))*tan(degtorad*(*pdec)) );
         *pdec = *pdec + 0.5*an*cos(degtorad*(*pra));
         
         printf("Pra, Pdec: %f %f\n",*pra,*pdec);
         
      }//end vertex loop

  return;
}
//////////////////////////////////////////////////////////////////////
// TEST: centroid
//////////////////////////////////////////////////////////////////////
int mainkl(int argc,char ** argv)
{
  float64 sumx=0, sumy=0, sum=0, I0=0;
  int sumxx=0, sumyy=0, summ=0;
  int i, j, xlen=0, ylen=0, xstart=0, ystart=0, ps=768000;
  char * nameFile1 = "TMO420001.nic";
  char * nameoutFile = "Test_C_out.grd";
  FILE  *inFile1, *outFile;
  short head1[3], head2[3];
  size_t itemsread, newlen;


  // Read in both files
  inFile1 = fopen(nameFile1,"rb");
  itemsread = fread(&head1, sizeof(short),3,inFile1);  
  itemsread = fread(&adata2, sizeof(short),ps,inFile1);
  fclose(inFile1);
  
  for(i=0; i<600; i++)
  {
    for(j=0; j<1280; j++)
    {
      adata4[i][j] = 0;
    }
  }
  
  for(i=80; i<145; i++)
  {
    for(j=675; j<740; j++)
    {
      adata4[i][j] = adata2[i][j];
      if(adata4[i][j]<0)
      {
        adata4[i][j] = adata4[i][j] + 65536;
      }
    }
  }
  
  I0 = 25000;
  
  for(i=80; i<145; i++)
  {
    for(j=675; j<740; j++)
    {
      adata4[i][j] = adata4[i][j] - I0;
      if(adata4[i][j]<0)
      {
        adata4[i][j] = 0;
      }

    }
  }

  xstart = 675;
  xlen = 740;
  ystart = 80; 
  ylen = 145;
  I0 = 0;  

     ////  Centroid algorithm  \\\\
  
  //xlen, ylen are 1 greater than there actual value{1->(len+1)}
  for(i=ystart; i<ylen; i++)
  {
    for(j=xstart; j<xlen; j++)
    {
      sum = sum + ( adata4[i][j] - I0 );
    }
  }
    
  for(i=xstart; i<xlen; i++)
  {
    for(j=ystart; j<ylen; j++)
    {
      sumx = sumx + ( i*(adata4[j][i] - I0) );
    }
  }

  for(i=ystart; i<ylen; i++)
  {
    for(j=xstart; j<xlen; j++)
    {
      sumy = sumy + ( i*(adata4[i][j] - I0) );
    }
  }
  

  printf("SumX, SumY, Sum: %f %f %f\n",sumx,sumy,sum);
  
  sumx = sumx/sum;
  sumy = sumy/sum;

  printf("Centroid: [%f,%f]\n",sumx, sumy);
  
  //adata4[119][709] = 0;

  /*
  printf("Writing to file: 'Test_C_out.grd'\n");
  
  //open file
  outFile = fopen(nameoutFile,"w");
      
  //write header
  fprintf(outFile,"DSAA\n1280 600\n0 1280\n0 600\n0 66000\n");

  //write array
  for (i=0; i<600;)
  {
      for (j=0; j<1280;)
      {
        fprintf(outFile,"%d ",adata4[i][j]);              
        j++;
      }
    i++;
  }
  fclose(outFile);
  */
  
  

  return 0;
}
//////////////////////////////////////////////////////////////////////
// TEST: find best ra0 adjustment for night 42xxxx
//////////////////////////////////////////////////////////////////////
int mainas(int argc,char ** argv)
{
  int i, j, k, l, h[8], ps=768000, *htemp;
  float64 aq, bq, cq, rq, thetax, thetaxsq, thetay, cx, cy, cz, crsq, cr, phi, rlat;
  float64 gamma=0.0, beta, alfa=0.50, camlat=34.9471, cra, cdec;
  float64 degperminute=0.2506844731,degtorad=0.0174533, radtodeg=57.29578;
  float64 x0=630.6, y0=1291.38, r0=1194.5, rccd, ra0;
  float64 response=0.0, x[8], y[8], r[8], d[8], noderesponse, *pra, *pdec;
  float64 ratotal=0, dectotal=0, rafinal=0, decfinal=0;
  float64 am=1.2812323, an=0.556753, t;
  float64 raoff=0, decoff=0, tempr, tempd;

          
       // adjustment for the night 42xxx
       ra0=279.425+0.97;//**ad hoc**
       
       //preession parameters from AQ4, p670
       t=0.017215;//time from J2000 in hundreds of y, for midnight 9/20/01
       am = am*t;
       an = an*t;
                                            
       x[0] = 124;
       y[0] = 205;
       x[1] = 488;
       y[1] = 101;
       x[2] = 509;
       y[2] = 101;
       x[3] = 1121;
       y[3] = 196;
       x[4] = 1070;
       y[4] = 177;
       x[5] = 645;
       y[5] = 93;
       x[6] = 393;
       y[6] = 117;
       x[7] = 783;
       y[7] = 113;
       
       h[0] = 180;
       h[1] = 260;
       h[2] = 325;
       h[3] = 325;
       h[4] = 357;
       h[5] = 386;
       h[6] = 443;
       h[7] = 443;
       
       r[0] = 326.046;
       d[0] = 9.875;
       r[1] = 345.943;
       d[1] = 28.082;
       r[2] = 2.097;
       d[2] = 29.091;
       r[3] = 2.294;
       d[3] = 59.15;
       r[4] = 10.127;
       d[4] = 56.539;
       r[5] = 17.433;
       d[5] = 35.621;
       r[6] = 31.793;
       d[6] = 23.463;
       r[7] = 30.975;
       d[7] = 42.330;
       
       //loop for each vertex of a pixel
       for (k=0; k<1; k++)
       {         
         rccd = sqrt((x[k] - x0)*(x[k] - x0) + (y[k] - y0)*(y[k] - y0));
         thetax = radtodeg*asin( (x[k] - x0) / rccd );
         thetax = thetax / 0.999744;//atmosphere adjustment
                        
         //quadratic coefficiants for the formula: rq=aqTY^2+bqTY+cq
         //from July 11, 2001 memo(ignoring 3rd and 4th order terms)
         /*
         thetaxsq = thetax*thetax;         
         aq = 0.1164 - ( pow(10,-5)*thetaxsq );
         bq = ( 3.6*pow(10,-5)*thetaxsq ) - 18.0266;
         cq = ( 0.000936*thetaxsq ) - ( 1.044*pow(10,-7)*thetaxsq*thetaxsq );
         */
         rq = (rccd - r0)/1.013;//y-plate scale adjustment  
         //thetay = ( -bq - sqrt( (bq*bq)-(4*aq*(cq-rq))) ) / ( 2*aq );
         thetay = rq / 18.0226; //note: +18.0226 to flip thetay/cy sign relative to Andy's program
                                                            
         cx = -sin(thetax*degtorad);
         cy = sin(thetay*degtorad);         
         crsq = (cx*cx)+(cy*cy);
         if(crsq>0)
         {
           phi = atan2(cy,cx);
         }else
         {
           phi = 0;
         }
         cz = sqrt(1-(crsq));
         rlat = radtodeg*asin(cz);
         phi = radtodeg*phi; 
                 
         beta = camlat-90;
         gamma = (-h[k]*degperminute)-ra0; //(frame#+1)*degperminute  
         pra = &phi;
         pdec = &rlat;          
         
         rotate(alfa, beta, gamma, pra, pdec);  
                 
         //precession adjustment
         *pra = *pra + 0.5*( am+an*sin(degtorad*(*pra))*tan(degtorad*(*pdec)) );
         *pdec = *pdec + 0.5*an*cos(degtorad*(*pra));
         
         tempr=*pra;
         tempd=*pdec;
         
         tempr = r[k] - tempr;
         if(tempr<0) tempr=tempr*(-1.0);
         tempd =d[k] - tempd;
         if(tempd<0) tempd=tempd*(-1.0);
         raoff = raoff + tempr;
         decoff = decoff + tempd;
         
         
         //printf("Ra, Dec: %f %f\n",*pra, *pdec);
         
       }
       
       //raoff = raoff/8;
       //decoff = decoff/8;
       
       printf("RAoff: %f\n",raoff);
       printf("DECoff: %f\n",decoff);
         

  return;
}











//////////////////////////////////////////////////////////////////
//	Main test Program
//////////////////////////////////////////////////////////////////
int mainop(int argc,char ** argv)
{
  char * nameoutFile = "Enif_T_pict.grd";
  char nameFile[] = "TMO420000.nic";
  FILE  *inFile, *outFile;
  short head[3];
  size_t itemsread, newlen;
  float64 ped=0, rup=0, rdown=360, dup=0, ddown=90;
  int pedcnt=0;
  int xarg, yarg, one=1, imlen, flen, fstart;
  int i, j, k, l, f, ps=768000, *htemp;
  float64 aq, bq, cq, rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq, cr, phi, rlat;
  float64 gamma=0.0, beta, alfa=0.50, camlat=34.9471, cra, cdec;
  float64 degperminute=0.2506844731,degtorad=0.0174533, radtodeg=57.29578;
  float64 x0=630.6, y0=1291.38, r0=1194.5, rccd, ra0=279.425;
  float64 response=0.0, x[4], y[4], noderesponse, *pra, *pdec;
  float64 *rtemp, totalarea=0, nodearea;
  SpatialIndex *si;
  SpatialVector *v[4], *tempv;
  SpatialConvex *pixelconvex;
  SpatialMarkup *markup;
  SpatialIDVec *partialvec, *fullvec, *imagevec;
  HitList *hlist;
  ResponseList *rlist;
  IDTYPE *partialnode, *fullnode, *imagenode;  
  
  printf("Starting...\n");
  
  // adjustment for the night 42xxx
  ra0=279.425 + 0.97; //**ad hoc**

  for(i=0; i<1600; i++)
  {
    for(j=0; j<400; j++)
    {
      image[i][j] = 0;
      imagecnt[i][j] = 0;
    }
  }
  
  //initialize for indexing
  si = spatialIndexNew( 12 , 2 );
  markup = spatialMarkupNew(si);    
  partialvec = (SpatialIDVec*) malloc( sizeof(SpatialIDVec));
  fullvec = (SpatialIDVec*) malloc( sizeof(SpatialIDVec));
  spatialIDVecInit(partialvec);
  spatialIDVecInit(fullvec);
  
  //Initialize the lists
  imagevec = (SpatialIDVec*) malloc( sizeof(SpatialIDVec));  
  hlist = (HitList*) malloc( sizeof(HitList));  
  rlist = (ResponseList*) malloc( sizeof(ResponseList));
  spatialIDVecInit(imagevec);
  hitListInit(hlist);
  responseListInit(rlist);
 
  printf("Looping...\n");
  
for (f=180; f<181; f++)
{    
  sprintf(&nameFile[5],"%04d.nic",f);
  //printf("inFile Name: %s\n",nameFile);
    
  for(i=0; i<600; i++)
  {
    for(j=0; j<1280; j++)
    {
      adata4[i][j] = 0;
      adata2[i][j] = 0;
    }
  }
  
  // Read in both file
  inFile = fopen(&nameFile,"rb");
  itemsread = fread(&head, sizeof(short),3,inFile);  
  itemsread = fread(&adata2, sizeof(short),ps,inFile);
  fclose(inFile);

  for(i=0; i<600; i++)
  {
    for(j=0; j<1280; j++)
    {
      adata4[i][j] = adata2[i][j];
      if(adata4[i][j]<0)
      {
        adata4[i][j] = adata4[i][j] + 65536;
      }
    }
  }

  //pedestal subtraction from both files
  for(i=2; i<600; i++)
  {
    for(j=1271; j<1280; j++)
    {
      ped = ped + adata4[i][j];
      pedcnt++;
    }
    for(j=4; j<13; j++)
    {
      ped = ped + adata4[i][j];
      pedcnt++;
    }
  }    
  ped = ped/pedcnt;     
  for(i=0; i<600; i++)
  {
    for(j=0; j<1280; j++)
    {
      adata4[i][j] = adata4[i][j] - ped;
    }
  }
  ped=0;
  pedcnt=0;
  
  //loop for each pixel
  for (i=0; i<600; i++)
  {
    for (j=0; j<1280; j++)
    {                            
     rccd = sqrt((j-x0)*(j-x0) + (i-y0)*(i-y0));
     dr = abs(rccd - r0);
     //less than 1.5 deg x 18.0226 pix/deg
     if( dr < 27.0339 )
     {
      if( abs((j-x0)/(i-y0)) < 0.57735 )// less than tan 30
      { 
       response = adata4[i][j];        
                                                    
       x[0] = j - 0.5;
       y[0] = i - 0.5;
       x[1] = j - 0.5;
       y[1] = i + 0.5;
       x[2] = j + 0.5;
       y[2] = i + 0.5;
       x[3] = j + 0.5;
       y[3] = i - 0.5;
       
       //loop for each vertex of a pixel
       for (k=0; k<4; k++)
       {         
         rccd = sqrt((x[k] - x0)*(x[k] - x0) + (y[k] - y0)*(y[k] - y0));
         thetax = radtodeg*asin( (x[k] - x0) / rccd );
         thetax = thetax / 0.999744;//atmosphere adjustment
                        
         //quadratic coefficiants for the formula: rq=aqTY^2+bqTY+cq
         //from July 11, 2001 memo(ignoring 3rd and 4th order terms)
         /*
         thetaxsq = thetax*thetax;         
         aq = 0.1164 - ( pow(10,-5)*thetaxsq );
         bq = ( 3.6*pow(10,-5)*thetaxsq ) - 18.0266;
         cq = ( 0.000936*thetaxsq ) - ( 1.044*pow(10,-7)*thetaxsq*thetaxsq );
         */
         rq = (rccd - r0)/1.013;//y-plate scale adjustment  
         //thetay = ( -bq - sqrt( (bq*bq)-(4*aq*(cq-rq))) ) / ( 2*aq );
         thetay = rq / 18.0226; //note: +18.0226 to flip thetay/cy sign relative to Andy's program
                                                            
         cx = -sin(thetax*degtorad);
         cy = sin(thetay*degtorad);         
         crsq = (cx*cx)+(cy*cy);
         if(crsq>0)
         {
           phi = atan2(cy,cx);
         }else
         {
           phi = 0;
         }
         cz = sqrt(1-(crsq));
         rlat = radtodeg*asin(cz);
         phi = radtodeg*phi; 
                 
         beta = camlat-90;
         gamma = (-f*degperminute)-ra0; //(frame#+1)*degperminute  
         pra = &phi;
         pdec = &rlat;          
         
         rotate(alfa, beta, gamma, pra, pdec); 
         
         //if(rup<*pra) rup=*pra;
         //if(rdown>*pra) rdown=*pra;
         //if(dup<*pdec) dup=*pdec;
         //if(ddown>*pdec) ddown=*pdec;
                 
         v[k] = spatialVectorNewRaDec(*pra, *pdec);
      }//end vertex loop
           
     pixelconvex = spatialConvexNewRectangle(v[0], v[1], v[2], v[3]);
     
     //free vector memory
     spatialVectorDelete(v[0]);
     spatialVectorDelete(v[1]);
     spatialVectorDelete(v[2]);
     spatialVectorDelete(v[3]);
          
     //intersect with index
     spatialConvexIntersect( pixelconvex, si, markup, partialvec, fullvec );
     
     //test center of partial nodes to see if they are within pixelconvex
     partialnode = partialvec->vector_;     
     for (k=0; k<partialvec->length_; k++, partialnode++)
     {        
        WspatialTriangleCenter( si, partialnode, pra, pdec );        
        tempv = spatialVectorNewRaDec(*pra, *pdec);       
        triangle[k] = spatialConvexTestVertex(pixelconvex,tempv);
        spatialVectorDelete(tempv);
        // triangle[k] = 1 or 0 for in or out of the pixel
        if (triangle[k] = 1) spatialIDVecAppend(fullvec, *partialnode);                
     }
     spatialConvexDelete(pixelconvex);
     
     //give each triangle an equal portion of the response
     fullnode = fullvec->vector_; 
     flen = fullvec->length_; 
     noderesponse = response/flen;
     
/*     for (k=0; k<flen; k++, fullnode++)
     {
       nodearea = spatialIndexAreaID(si, *fullnode);
       totalarea = totalarea + nodearea;
     }   
     fullnode = fullvec->vector_;
*/     
     imlen = imagevec->length_;
     fstart = 0;
     if(imlen==0) 
     {
//       nodearea = spatialIndexAreaID(si, *fullnode);
//       noderesponse = response*(nodearea/totalarea);
       
       spatialIDVecAppend(imagevec, *fullnode);
       responseListAppend(rlist, noderesponse);
       hitListAppend(hlist, one);
       fullnode++;
       fstart++;
       imlen = imagevec->length_;//new length after appendation
     }           
     for (k=fstart; k<flen; k++, fullnode++)
     {
       imagenode = imagevec->vector_;
       htemp = hlist->hits_;
       rtemp = rlist->response_;       
       for (l=0; l<imlen; l++, imagenode++, htemp++, rtemp++)
       {
         if(*fullnode == *imagenode)
         {
//           nodearea = spatialIndexAreaID(si, *fullnode);
//           noderesponse = response*(nodearea/totalarea);
           *rtemp = *rtemp + noderesponse;
           *htemp = *htemp + 1;
           break;
         }
         if(l == (imlen - 1) )
         {
//          nodearea = spatialIndexAreaID(si, *fullnode);
//           noderesponse = response*(nodearea/totalarea);
           spatialIDVecAppend(imagevec, *fullnode);
           responseListAppend(rlist, noderesponse);
           hitListAppend(hlist, one);
           break;
         }
       }
     }    
     totalarea = 0;
     
     // clear IDvecs    
     spatialIDVecClear(partialvec, true);
     spatialIDVecClear(fullvec, true);

      
     }//end pixel's if #1
    }//end pixel's if #2
   }//end pixel loop(j)   
   printf("%d: loop Index(i): %d\n",f,i);   
  }//end pixel loop(i)
}//end f frame loop

  //free IDvec & markup memory
  free(partialvec);
  free(fullvec);
  free(markup);
 
  //approximate center of FOV
  cra=325;
  cdec=35;
  
 
  //printf("Length, Capacity: %d %d\n",imagevec->length_,imagevec->capacity_);    
  printf("Imaging...\n");
  
  // loop to put all triangles with averaged response into data array
  imagenode = imagevec->vector_;
  rtemp = rlist->response_;  
  htemp = hlist->hits_;
  for(k=0; k<imagevec->length_; k++, imagenode++, htemp++, rtemp++)
  {   
    noderesponse = (*rtemp) / (*htemp);
    WspatialTriangleCenter( si, imagenode, pra, pdec );
    *pra = (*pra - cra)*cos(*pdec*degtorad) + cra;

    xarg =  ( (*pra) - cra )*20 +199;//x-coordinate
    yarg =  ( (*pdec) - cdec )*20 + 799;//y-coordinate        
    //make sure xarg/yarg are within array bounds
    if (xarg<0)
    {
      xarg=0;
    }else if (xarg>=400)
    {
      xarg=399;
    }
    if (yarg<0)
    {
      yarg=0;
    }else if (yarg>=1600)
    {
      yarg=1599;
    }
      
    //image[y][x] because C fills columns before rows
    image[yarg][xarg] = image[yarg][xarg] + noderesponse;
    imagecnt[yarg][xarg] = imagecnt[yarg][xarg] + 1;
  }//end data array loop

/*  
  for(i=0; i<300; i++)
  {
    for(j=0; j<300; j++)
    {
    if (imagecnt[i][j] > 0)  image[i][j] = image[i][j]/imagecnt[i][j];      
    }
  }
*/

  printf("Writing 'Enif_T_pict.grd'\n");
  
  //loop for writing the data array to a file
  outFile = fopen(nameoutFile,"w");
      
  //write header
  fprintf(outFile,"DSAA\n1600 400\n315 335\n-5 75\n0 6000\n");

  //write array
  for (i=0; i<1600;)
  {
      for (j=0; j<400;)
      {       
       fprintf(outFile,"%8.6f ",image[i][j]);
                  
       j++;
      }
    i++;
  }
  fclose(outFile);
  
  
  //printf("%f < RA < %f\n",rdown,rup);
  //printf("%f < DEC < %f\n",ddown,dup);

  printf("Finished!!!\n");

  return 0;
}







//test program to get the indices of the two example nodes from lookup.c
int mainwes(int argc, const char *argv[])
{
  short index[INDEX_SIZE];
  int i,j;
  
  unsigned int ids[2] = {2486622255, 16020};//from lookup.c
  
  for(i=0; i<2; i++)
  {
    printf("\nid: %u\n",ids[i]);
    printf("Indices: \n");
    
    getIndices(ids[i], index);
    
    for(j=0; j<16; j++)
    {
      printf("%d ",index[j]);
    }
  }
  
  return 0;
}


int mainjaj(int argc,char ** argv) 
{
  SpatialIndex *si;
  SpatialVector *v[5];
  SpatialConvex *pixelconvex;
  SpatialMarkup *markup;
  SpatialIDVec *partialvec, *fullvec;
  float64 ra[5], dec[5];
  int k, i, j;
  time_t t1, t2, t=0;
  
  //initialize for indexing
  si = spatialIndexNew( 12 , 2 );
  markup = spatialMarkupNew(si);    
  partialvec = (SpatialIDVec*) malloc( sizeof(SpatialIDVec));
  fullvec = (SpatialIDVec*) malloc( sizeof(SpatialIDVec));
  spatialIDVecInit(partialvec);
  spatialIDVecInit(fullvec);
  
  ra[1]=328.95;
  ra[2]=328.90;
  ra[3]=328.88;
  ra[4]=328.93;
  dec[1]=10.67;
  dec[2]=10.66;
  dec[3]=10.70;
  dec[4]=10.71;

  for(k=1; k<5; k++)
  {
    v[k] = spatialVectorNewRaDec(ra[k], dec[k]);
  }
  pixelconvex = spatialConvexNewRectangle(v[1], v[2], v[3], v[4]);
   
  for(j=0; j<20; j++)
  { 
   t1=clock();
   for(i=0; i<1000; i++)
   {      
    spatialConvexIntersect( pixelconvex, si, markup, partialvec, fullvec );
    spatialIDVecClear(fullvec, true);
    spatialIDVecClear(partialvec, true);
   }
   t2=clock();
   t=t2-t1;
   printf("%d Time: %f\n",j,(double)t/(double)CLOCKS_PER_SEC);
   fflush(stdout);
  }

  return 0;
}