#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "obs.h"

/*.......................................................................
 * Determine the range of UV radii, |U|, and |V| distances, and
 * amplitudes, of visibilities in the current IF and within the UV
 * radius range uvmin -> uvmax (inclusive).
 *
 * Input:
 *  ob     Observation *  The observation to be parameterised.
 *  doall          int    Find the range for all IFs if true. Otherwise
 *                        find the range from just the current stream IF.
 *  uvmin        float    The UV radius (wavelengths) below which to ignore
 *                        data.
 *  uvmax        float    The UV radius (wavelengths) beyond which to ignore
 *                        data. If the largest of uvmin and uvmax is <= 0.0f
 *                        then the returned limits pertain to the whole data
 *                        set.
 * Output:
 *  return    UVrange *   Pointer to a static internal container, recording
 *                        the maximum UV radius, |U|,|V|, and amplitude of
 *                        Visibilities within UV radius uvmin to uvmax.
 *                        On error, NULL is returned.
 */
UVrange *uvrange(Observation *ob, int doall, float uvmin, float uvmax)
{
  static UVrange uvr; /* The container of statistics to be returned */
  int ut;        /* The ut being looked at */
  int base;      /* The baseline being looked at */
  int docut;     /* If true, apply uvcut UV radius cutoff */
  int cif;       /* The index of the IF being processed */
  int bif,eif;   /* The indexes of the first and last index to be processed */
  float minrsq;  /* Min UV radius squared of a visibility yet encounterred */
  float maxrsq;  /* Max UV radius squared of a visibility yet encounterred */
  int first = 1; /* True until the first in-range visibility is encounterred */
  int old_if;    /* The index of IF to be restored on return */
/*
 * Check the observation state.
 */
  if(!ob_ready(ob, doall ? OB_SELECT : OB_GETIF, "uvrange"))
    return NULL;
/*
 * Record the index of the current IF.
 */
  old_if = get_cif_state(ob);
/*
 * Enforce positivity.
 */
  if(uvmin < 0.0f)
    uvmin = 0.0f;
  if(uvmax < 0.0f)
    uvmax = 0.0f;
/*
 * Arrange that uvmin <= uvmax.
 */
  if(uvmin > uvmax) {float ftmp = uvmin; uvmin = uvmax; uvmax = ftmp;};
/*
 * Record whether a UV cutoff range was specified and precalculate the
 * radius squared of uvmin and uvmax to remove the requirement for a
 * square root in sqrt(U^2+V^2) when comparing real UV radii to the
 * range limits.
 */
  docut = uvmax > 0.0f;
  uvmin *= uvmin;
  uvmax *= uvmax;
/*
 * Initialize the container.
 */
  minrsq = maxrsq = 0.0f;
  uvr.uvrmin = uvr.uvrmax = 0.0f;
  uvr.umin = uvr.umax = 0.0f;
  uvr.vmin = uvr.vmax = 0.0f;
  uvr.ampmin = uvr.ampmax = 0.0f;
/*
 * Set the range of IFs to be processed.
 */
  if(doall) {
    bif = 0;
    eif = ob->nif-1;
  } else {
    bif = eif = ob->stream.cif;
  };
/*
 * Loop through all IFs, ending with the default IF.
 */
  for(cif=bif; (cif=nextIF(ob, cif, 1, 1))>=0 && cif<=eif; cif++) {
    float uvscale;  /* UV coordinate scale factor */
    Subarray *sub;  /* The descriptor of the sub-array being processed */
    int isub;       /* The index of sub in ob->sub[] */
/*
 * Get the next IF.
 */
    if(getIF(ob, cif))
      return NULL;
/*
 * Get the uvscale factor for the new IF.
 */
    uvscale = ob->stream.uvscale;
/*
 * Look at all integrations of all sub-arrays.
 */
    sub = ob->sub;
    for(isub=0; isub<ob->nsub; isub++, sub++) {
      Integration *integ = sub->integ;
/*
 * Update maxrsq and maxamp for each visibility that exceeds it.
 */
      for(ut=0; ut<sub->ntime; ut++,integ++) {
	Visibility *vis = integ->vis;
	for(base=0; base<sub->nbase; base++,vis++) {
/*
 * Ignore flagged visibilities.
 */
	  if(!vis->bad) {
/*
 * Get the U and V coordinates of the visbility.
 */
	    float uu = fabs(vis->u * uvscale);
	    float vv = fabs(vis->v * uvscale);
/*
 * Get the UV radius squared.
 */
	    float newrsq = uu*uu + vv*vv;
/*
 * Ignore all visibilities except those between uvmin and uvmax.
 */
	    if(!docut || (newrsq >= uvmin && newrsq <= uvmax)) {
/*
 * Initialize the ranges with the first visibility that lies in range.
 */
	      if(first) {
		first = 0;
		uvr.ampmin = uvr.ampmax = vis->amp;
		minrsq = maxrsq = newrsq;
		uvr.umin = uvr.umax = uu;
		uvr.vmin = uvr.vmax = vv;
	      } else {
/*
 * Update the amplitude range.
 */
		if(vis->amp > uvr.ampmax)
		  uvr.ampmax = vis->amp;
		if(vis->amp < uvr.ampmin)
		  uvr.ampmin = vis->amp;
/*
 * Update the min/max UV radius.
 */
		if(newrsq < minrsq)
		  minrsq = newrsq;
		else if(newrsq > maxrsq)
		  maxrsq = newrsq;
/*
 * Update the min/max U distance.
 */
		if(uu < uvr.umin)
		  uvr.umin = uu;
		else if(uu > uvr.umax)
		  uvr.umax = uu;
/*
 * Update the min/max V distance.
 */
		if(vv < uvr.vmin)
		  uvr.vmin = vv;
		else if(vv > uvr.vmax)
		  uvr.vmax = vv;
	      };
	    };
	  };
	};
      };
    };
  };
/*
 * Convert the min/max UV radii squared to their equivalent radii.
 */
  uvr.uvrmin = sqrt(minrsq);
  uvr.uvrmax = sqrt(maxrsq);
/*
 * Restore the start IF.
 */
  if(set_cif_state(ob, old_if))
    return NULL;
/*
 * Return the container that holds the results.
 */
  return &uvr;
}