#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "logio.h"
#include "vlbconst.h"
#include "model.h"
#include "besj.h"

#define CTSIZ 4096                 /* Size of lookup table */

static float costab[CTSIZ+2];      /* Cosine lookup table */
static float *cos_tab=&costab[1];  /* Allow elements -1 to CTSIZ */
static int ready=0;                /* True after table is initialized */
static const int soffset=CTSIZ+(CTSIZ/4); /* Index to cos(360+90 degrees) */

/*.......................................................................
 * Compute the visibility amplitude and phase at a given U,V coordinate
 * corresponding to a list of model components.
 *
 * Input:
 *   mod  Model *  The input model.
 *   uu   float    The U coordinate in the UV plane (wavelengths).
 *   vv   float    The V coordinate in the UV plane (wavelengths).
 * Output:
 *   amp  float *  The corresponding amplitude.
 *   phs  float *  The corresponding phase (radians).
 */
void modvis(Model *mod, float uu, float vv, float *amp, float *phs)
{
  float imvis=0.0;  /* The accumulated imaginary component of the visibility */
  float revis=0.0;  /* The accumulated real component of the visibility */
  Modcmp *cmp;      /* The model component being used */
  float cmpphs;     /* Component phase */
  float cmpamp;     /* Component amplitude */
  float ftmp;       /* Work variables */
  float *cos_ptr;   /* Pointer into cosine lookup table */
  float err_indx;   /* Holds quantisation error in a lookup table index */
  int isign;        /* Sign of an angle */
  int cos_indx;     /* Index into cosine lookup table */
/*
 * Initialise cosine lookup table on first entry to this function.
 */
  if(!ready) {
    ready=1;
    for(cos_indx=-1; cos_indx<CTSIZ+1; cos_indx++)
      cos_tab[cos_indx]=cos(twopi*cos_indx/CTSIZ);
  };
/*
 * Handle models which are completely comprised of delta-components
 * in a separate optimized loop.
 */
  if(mod->isdelt) {
    for(cmp=mod->head; cmp != NULL; cmp=cmp->next) {
      cmpamp = cmp->flux;
      cmpphs = (uu*cmp->x + vv*cmp->y); /* Component phase divided by twopi */
/*
 * Use interpolation into the cosine lookup table to get cos(twopi*cmpphs).
 */
      isign = (cmpphs<0.0f)?-1:1;     /* Record sign of angle */
      ftmp = cmpphs*(isign*CTSIZ);    /* Decimal index into lookup table */
      cos_indx = (int) ftmp;          /* Integer index into lookup table */
      err_indx = ftmp - cos_indx;     /* Quantisation error in cos_indx */
      cos_indx %= CTSIZ;              /* Fold index into twopi range */
      cos_ptr = &cos_tab[cos_indx];   /* Table element of truncated index */
/*
 * Finally, compute revis += cmpamp*cos(twopi*cmpphs).
 */
      revis += cmpamp * (*cos_ptr + err_indx*(cos_ptr[1] - *cos_ptr));
/*
 * Now use interpolation to determine sine(twopi*cmpphs).
 */
      cos_indx = (soffset - isign*cos_indx) % CTSIZ; /* 360+(90-angle) */
      cos_ptr  = &cos_tab[cos_indx];
/*
 * And compute imvis += cmpamp*sin(twopi*cmpphs).
 */
      imvis += cmpamp * (*cos_ptr + err_indx*(cos_ptr[-isign] - *cos_ptr));
    };
  } else {
/*
 * Models with mixed non-delta and delta model components.
 */
    for(cmp=mod->head; cmp != NULL; cmp=cmp->next) {
/*
 * All the components are even functions so the phase is just the
 * fourier component phase at the centroid of the component.
 */
      cmpphs = (float) twopi * (uu*cmp->x + vv*cmp->y);
/*
 * The component amplitude is different for each component type.
 * Delta components are the easiest and commonest.
 */
      if(cmp->type == M_DELT) {
	cmpamp = cmp->flux;
      }
/*
 * The other types are a little more complicated!
 */
      else {
/*
 * Pre-compute parameters.
 */
	double sinphi = sin(cmp->phi);
	double cosphi = cos(cmp->phi);
	double tmpa = (vv*cosphi+uu*sinphi);
	double tmpb = (cmp->ratio*(uu*cosphi-vv*sinphi));
	double tmpc = pi * cmp->major * sqrt(tmpa*tmpa + tmpb*tmpb);
/*
 * Limit tmpc to sensible values to prevent underflow,overflow and
 * divide-by-zero errors.
 */
	if(tmpc < 1.0e-9)
	  tmpc = 1.0e-9;
/*
 * See the "Introduction to Caltech VLBI programs" for details of
 * the type-specific calculations below.
 */
	switch(cmp->type) {
	case M_GAUS:
	  cmpamp = cmp->flux * (tmpc<12.0 ? exp(-0.3606737602 * tmpc*tmpc):0.0);
	  break;
	case M_DISK:
	  cmpamp = 2.0f * cmp->flux * c_besj1(tmpc)/tmpc;
	  break;
	case M_ELLI:
	  cmpamp = 3.0f * cmp->flux * (sin(tmpc)-tmpc*cos(tmpc))/(tmpc*tmpc*tmpc);
	  break;
	case M_RING:
	  cmpamp = cmp->flux * c_besj0(tmpc);
	  break;
	case M_RECT:
	  tmpa = pi * cmp->major * (uu*sinphi+vv*cosphi);
	  cmpamp = cmp->flux * (fabs(tmpa) > 0.001 ? sin(tmpa)/tmpa : 1.0);
	  break;
	case M_SZ:
	  cmpamp = cmp->flux * (tmpc < 50.0 ? exp(-tmpc) : 0.0) / tmpc;
	  break;
	default:
	  lprintf(stderr, "Ignoring unknown model component type: %d\n",
		  cmp->type);
	  continue;
	};
      };
/*
 * Accuumulate real and imaginary parts.
 */
      revis += cmpamp*cos(cmpphs);
      imvis += cmpamp*sin(cmpphs);
    };
  };
/*
 * Convert the accumulated complex visibility to amplitude and phase.
 */
  *amp = sqrt(imvis*imvis + revis*revis);
  *phs = (imvis==0.0f && revis==0.0f) ? 0.0f : atan2(imvis,revis);
  return;
}

/*.......................................................................
 * Compute the visibility amplitude and phase at a given U,V coordinate
 * corresponding to a single model component.
 *
 * Input:
 *   cmp Modcmp *  The descriptor of the model component.
 *   uu   float    The U coordinate in the UV plane (wavelengths).
 *   vv   float    The V coordinate in the UV plane (wavelengths).
 * Output:
 *   amp  float *  The corresponding amplitude.
 *   phs  float *  The corresponding phase (radians).
 */
void cmpvis(Modcmp *cmp, float uu, float vv, float *amp, float *phs)
{
  float cmpamp;  /* The amplitude of the component */
/*
 * All the components are even functions so the phase is just the
 * fourier component phase at the centroid of the component.
 */
  float cmpphs = (float) twopi * (uu*cmp->x + vv*cmp->y);
/*
 * The component amplitude is different for each component type.
 * Delta components are the easiest and commonest.
 */
  if(cmp->type == M_DELT) {
    cmpamp = cmp->flux;
  }
/*
 * The other types are a little more complicated!
 */
  else {
/*
 * Pre-compute parameters.
 */
    double sinphi = sin(cmp->phi);
    double cosphi = cos(cmp->phi);
    double tmpa = (vv*cosphi+uu*sinphi);
    double tmpb = (cmp->ratio*(uu*cosphi-vv*sinphi));
    double tmpc = pi * cmp->major * sqrt(tmpa*tmpa + tmpb*tmpb);
/*
 * Limit tmpc to sensible values to prevent underflow,overflow and
 * divide-by-zero errors.
 */
    if(tmpc < 1.0e-9)
      tmpc = 1.0e-9;
/*
 * See the "Introduction to Caltech VLBI programs" for details of
 * the type-specific calculations below.
 */
    switch(cmp->type) {
    case M_GAUS:
      cmpamp = cmp->flux * (tmpc<12.0 ? exp(-0.3606737602 * tmpc*tmpc):0.0);
      break;
    case M_DISK:
      cmpamp = 2.0f * cmp->flux * c_besj1(tmpc)/tmpc;
      break;
    case M_ELLI:
      cmpamp = 3.0f * cmp->flux * (sin(tmpc)-tmpc*cos(tmpc))/(tmpc*tmpc*tmpc);
      break;
    case M_RING:
      cmpamp = cmp->flux * c_besj0(tmpc);
      break;
    case M_RECT:
      tmpa = pi * cmp->major * (uu*sinphi+vv*cosphi);
      cmpamp = cmp->flux * (fabs(tmpa) > 0.001 ? sin(tmpa)/tmpa : 1.0);
      break;
    case M_SZ:
      cmpamp = cmp->flux * (tmpc < 50.0 ? exp(-tmpc) : 0.0) / tmpc;
      break;
    default:
      lprintf(stderr, "Ignoring unknown model component type: %d\n", cmp->type);
      cmpamp = 0.0f;
      cmpphs = 0.0f;
    };
  };
/*
 * Assign the return values.
 */
  *amp = cmpamp;
  *phs = cmpphs;
  return;
}