/*****************************************************************************
******************************************************************************
**<AUTO>
**
** FILE:	area.c
**
**<HTML>
**	This file contains C functions to calculate survey areas of fields.
**</HTML>
**</AUTO>
**
**
** ENVIRONMENT:
**	ANSI C.
**
******************************************************************************
******************************************************************************
*/
#include <dervish.h>
#include "taCalibObj.h"

/*****************************************************************************
******************************************************************************
**<AUTO EXTRACT>
**
** ROUTINE:	taFieldArea
**
**<HTML>
**	Calculate the area associated with a field within the segment, both
**	for the OK_SCANLINE sample and the PRIMARY sample.
**</HTML>
**</AUTO>
******************************************************************************
******************************************************************************
*/
void
taFieldArea(const TRANS *trans,     /* IN: r' TRANS structure for the field */
	    const TA_GEOMETRY *geometry, /* IN: segment geometry */
	    double *segmentArea,    /* OUT: area of segment covered by the
				     *      field (degrees^2) */
	    double *primaryArea     /* OUT: primary area covered by the
				     *      field (degrees^2) */)
{
  double mu0, mu1, nu, nuMin, nuMax, tmp, buffer, nu0, nu1, dMu, muMid;
  double mu, nuStripe, nuHi, nuLambda, tanEta2, scale;
  int lambdaBound, nSteps, i;

  /* Verify proper geometry structure */
  shAssert((geometry->muLimits[0] < -999 && geometry->muLimits[1] > 999) ||
	   (geometry->muLimits[0] > -999 && geometry->muLimits[1] < 999));
  shAssert((geometry->lambdaLimits[0]<-999 && geometry->lambdaLimits[1]>999) ||
	   (geometry->lambdaLimits[0]>-999 && geometry->lambdaLimits[1]<999));
  shAssert((geometry->etaLimits[0] < -999 && geometry->etaLimits[1] > 999) ||
	   (geometry->etaLimits[0] > -999 && geometry->etaLimits[1] < 999));
  if (geometry->lambdaLimits[0] > -999)
    shAssert(geometry->etaLimits[0] > -999);
  if (geometry->etaLimits[0] > -999)
    shAssert(fabs(geometry->etaLimits[1] - geometry->etaLimits[0]
		  - at_stripeSeparation) < 0.0001);

  /* Fetch beginning and ending mu of the field (at center column) */
  atTransApply(trans, 'r', geometry->nOverlap / 2.,
	       0., geometry->nCols / 2., 0., NULL, NULL, &mu0, NULL, &nu,
	       NULL);
  atTransApply(trans, 'r', geometry->nRowsPerField + geometry->nOverlap / 2.,
	       0., geometry->nCols / 2., 0., NULL, NULL, &mu1, NULL, &nu,
	       NULL);

  /* Limit scan area to bounded region in mu.
   * Properly handle southern stripes, which wrap around zero. */
  if (geometry->muLimits[0] > -999)
    {
      if (geometry->muLimits[0] > 275.)
	{
	  if (mu0 < 275.) mu0 += 360.;
	  if (mu1 < 275.) mu1 += 360.;
	}
      mu0 = (geometry->muLimits[0] > mu0 ? geometry->muLimits[0] : mu0);
    }
  if (geometry->muLimits[1] < 999)
    {
      if (geometry->muLimits[1] > 275.)
	{
	  if (mu0 < 275.) mu0 += 360.;
	  if (mu1 < 275.) mu1 += 360.;
	}
      mu1 = (geometry->muLimits[1] < mu1 ? geometry->muLimits[1] : mu1);
    }
  if (mu0 > mu1) mu1 = mu0;
  
  /* Calculate scanline area. */
  scale = cos(nu * at_deg2Rad);
  *segmentArea = (mu1 - mu0) *
    (geometry->nuLimits[1] - geometry->nuLimits[0]) * scale;

  /* Convert mu to be zero at the pole. */
  mu0 -= 95.;
  mu1 -= 95.;

  /* If neither eta nor lambda limits apply, then stripe area is same
   * as scanline area.  Lambda limits are applied only when eta limits
   * are applied, and then only for skyVersion=0 (an assert above
   * guarantees that).  */
  if (geometry->etaLimits[0] < -999)
    {
      *primaryArea = *segmentArea;
      return;
    }

  /* Applying eta (and possibly lambda) cuts.
   * The northern stripe eta boundary can be expressed in (mu,nu) coords as:
   *     nu = tan(at_stripeSeparation/2.) * sin(mu)
   */
  tanEta2 = tan(at_stripeSeparation * at_deg2Rad / 2.);

  /* Determine min and max nu for the scanline. */
  nuMin = geometry->nuLimits[0];
  nuMax = geometry->nuLimits[1];
  /* Flip to the north, just to simplify stuff. */
  if (nuMin < 0 && fabs(nuMin) > fabs(nuMax))
    {
      tmp = nuMin;
      nuMin = -nuMax;
      nuMax = -tmp;
    }

  /* Does the stripe lambda limits cross the field? */
  if (geometry->lambdaLimits[0] > -999)
    {
      buffer = 25. * at_asec2Deg; /* Max change in mu along lambda limit
				   * across a stripe */
      if (mu1 < geometry->lambdaLimits[0] + 90 - buffer ||
	  mu0 > geometry->lambdaLimits[1] + 90 + buffer)
	{
	  /* Field completely beyond stripe lambda limits. */
	  *primaryArea = 0;
	  return;
	}
      lambdaBound = 0;
      if (geometry->lambdaLimits[0] + 90 > mu0 &&
	  geometry->lambdaLimits[0] + 90 - buffer < mu1)
	lambdaBound = -1;
      else if (geometry->lambdaLimits[1] + 90 + buffer > mu0
	       && geometry->lambdaLimits[1] + 90 < mu1)
	lambdaBound = 1;
    }
  else
    lambdaBound = 0;

  /* If the stripe lambda limits apply, then we must integrate. */
  if (lambdaBound != 0)
    {
      *primaryArea = 0;
      nSteps = 1000;
      dMu = (mu1 - mu0) / nSteps;
      for (i = 0; i < nSteps; i++)
	{
	  mu = mu0 + (i + 0.5) * dMu;
	  nuStripe = tanEta2 * sin(mu * at_deg2Rad) * at_rad2Deg;
	  if (nuStripe < nuMin) continue;
	  nuHi = (nuStripe < nuMax ? nuStripe : nuMax);
	  if (lambdaBound == -1)
	    {
	      if (mu >= geometry->lambdaLimits[0] + 90)
		*primaryArea += dMu * (nuHi - nuMin);
	      else
		{
		  if (mu < geometry->lambdaLimits[0] + 90 - buffer) continue;
		  nuLambda =
		    acos(cos((geometry->lambdaLimits[0] + 90) * at_deg2Rad) /
			 cos(mu * at_deg2Rad)) * at_rad2Deg;
		  if (nuLambda > nuHi) continue;
		  if (nuLambda < nuMin)
		    *primaryArea += dMu * (nuHi - nuMin) * scale;
		  else
		    *primaryArea += dMu * (nuHi - nuLambda) * scale;
		}
	    }
	  else
	    {
	      if (mu <= geometry->lambdaLimits[1] + 90)
		*primaryArea += dMu * (nuHi - nuMin);
	      else
		{
		  if (mu > geometry->lambdaLimits[1] + 90 + buffer) continue;
		  nuLambda =
		    acos(cos((geometry->lambdaLimits[1] + 90) * at_deg2Rad) /
			 cos(mu * at_deg2Rad)) * at_rad2Deg;
		  if (nuLambda > nuHi) continue;
		  if (nuLambda < nuMin)
		    *primaryArea += dMu * (nuHi - nuMin) * scale;
		  else
		    *primaryArea += dMu * (nuHi - nuLambda) * scale;
		}
	    }
	}
      return;
    }

  /* Stripe lambda limits don't matter. */
  /* Does the stripe eta boundary cross the field */
  nu0 = tanEta2 * sin(mu0 * at_deg2Rad) * at_rad2Deg;
  nu1 = tanEta2 * sin(mu1 * at_deg2Rad) * at_rad2Deg;
  if (nu0 >= nuMax && nu1 >= nuMax)
    {
      /* Entire scanline area in the stripe */
      *primaryArea = *segmentArea;
    }
  else if (nu0 <= nuMin && nu1 <= nuMin)
    {
      /* Entire scanline area outside the stripe */
      *primaryArea = 0;
    }
  else
    {
      /* Part of scanline in stripe area */
      if (nu0 >= nuMin && nu0 <= nuMax && nu1 >= nuMin && nu1 <= nuMax)
	{
	  /* Stripe crosses through entire field.  Keep same mu limits. */
	  *primaryArea = scale *
	    (-tanEta2 * (cos(mu1 * at_deg2Rad) - cos(mu0 * at_deg2Rad))
	     * at_rad2Deg * at_rad2Deg - nuMin * (mu1 - mu0));
	}
      else if (nu0 < nuMin || nu1 < nuMin)
	{
	  /* Enters or leaves below scanline */
	  muMid = asin(nuMin * at_deg2Rad / tanEta2) * at_rad2Deg;
	  if (muMid < mu0) muMid = 180. - muMid;
	  shAssert(muMid >= mu0 && muMid <= mu1);
	  if (nu0 < nuMin)
	    mu0 = muMid;
	  else
	    mu1 = muMid;
	  *primaryArea = scale *
	    (-tanEta2 * (cos(mu1 * at_deg2Rad) - cos(mu0 * at_deg2Rad))
	     * at_rad2Deg * at_rad2Deg - nuMin * (mu1 - mu0));
	}
      else
	{
	  /* Enters or leaves above scanline */
	  muMid = asin(nuMax * at_deg2Rad / tanEta2) * at_rad2Deg;
	  if (muMid < mu0) muMid = 180. - muMid;
	  shAssert(muMid >= mu0 && muMid <= mu1);
	  if (nu0 > nuMax)
	    mu0 = muMid;
	  else
	    mu1 = muMid;
	  *primaryArea = *segmentArea + scale *
	    ((-tanEta2 * (cos(mu1 * at_deg2Rad) - cos(mu0 * at_deg2Rad))
	      * at_rad2Deg * at_rad2Deg - nuMax * (mu1 - mu0)));
	}
    }
  return;
}