/*-
 * Written by H. Mitasova, L. Mitas, I. Kosinovsky, D. Gerdes Fall 1994
 * University of Illinois
 * US Army Construction Engineering Research Lab  
 * Copyright 1994, H. Mitasova (University of Illinois),
 * L. Mitas (University of Illinois),
 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)   
 *
 * modified by McCauley in August 1995
 * modified by Mitasova in August 1995  
 *
 */

#include <stdio.h>
#include <math.h>
#include <unistd.h>
#include <grass/gis.h>
#include <grass/bitmap.h>
#include <grass/interpf.h>

int IL_secpar_loop_2d (
    struct interp_params *params,
    int ngstc,			/* starting column */
    int nszc,			/* ending column */
    int k,			/* current row */
    struct BM *bitmask,
    double *gmin,
    double *gmax,
    double *c1min,
    double *c1max,
    double *c2min,
    double *c2max,	/* min,max interp.
							 * values */
    int cond1,
    int cond2		/* determine if particular values need to
				 * be computed */
)

/*
 * Computes slope, aspect and curvatures (depending on cond1, cond2) for
 * derivative arrays adx,...,adxy between columns ngstc and nszc.
 */
{
  double dnorm1, ro,		/* rad to deg conv */
   dx2=0, dy2=0, grad2=0,	/* gradient squared */
   slp=0, grad,			/* gradient */
   oor=0,			/* aspect  (orientation) */
   curn=0,			/* profile curvature */
   curh=0,			/* tangential curvature */
   curm=0,			/* mean curvature */
   temp,			/* temp  variable */
   dxy2;			/* temp variable   square of part diriv. */

  double gradmin;
  int i, got, bmask = 1;
  static int first_time_g = 1;

  ro = 57.295779;
  gradmin = 0.001;


  for (i = ngstc; i <= nszc; i++)
  {
    if (bitmask != NULL)
    {
      bmask = BM_get (bitmask, i, k);
    }
    got = 0;
    if (bmask == 1)
    {
      while ((got == 0) && (cond1))
      {
	dx2 = (double) (params->adx[i] * params->adx[i]);
	dy2 = (double) (params->ady[i] * params->ady[i]);
	grad2 = dx2 + dy2;
	grad = sqrt (grad2);
	/* slope in %        slp = 100. * grad; */
	/* slope in degrees */
	slp = ro * atan (grad);
	if (grad <= gradmin)
	{
	  oor = 0.;
	  got = 3;
	  if (cond2)
	  {
	    curn = 0.;
	    curh = 0.;
	    got = 3;
	    break;
	  }
	}
	if (got == 3)
	  break;

	/***********aspect from r.slope.aspect, with adx, ady computed
	            from interpol. function RST **************************/

	if (params->adx[i] == 0.)
	{
	  if (params->ady[i] > 0.)
	    oor = 90;
	  else
	    oor = 270;
	}
	else
	{
	  oor = ro * atan2 (params->ady[i], params->adx[i]);
	  if (oor <= 0.)
	    oor = 360. + oor;
	}

	got = 1;
      }				/* while */
      if ((got != 3) && (cond2))
      {

	dnorm1 = sqrt (grad2 + 1.);
	dxy2 = 2. * (double) (params->adxy[i] * params->adx[i] * params->ady[i]);


	curn = (double) (params->adxx[i] * dx2 + dxy2 + params->adyy[i] * dy2) / (grad2 * dnorm1 * dnorm1 * dnorm1);

	curh = (double) (params->adxx[i] * dy2 - dxy2 + params->adyy[i] * dx2) / (grad2 * dnorm1);

	temp = grad2 + 1.;
	curm = .5 * ((1. + dy2) * params->adxx[i] - dxy2 + (1. + dx2) * params->adyy[i]) / (temp * dnorm1);
      }
      if (first_time_g)
      {
	first_time_g = 0;
	*gmin = *gmax = slp;
	*c1min = *c1max = curn;
	*c2min = *c2max = curh;
      }
      *gmin = amin1 (*gmin, slp);
      *gmax = amax1 (*gmax, slp);
      *c1min = amin1 (*c1min, curn);
      *c1max = amax1 (*c1max, curn);
      *c2min = amin1 (*c2min, curh);
      *c2max = amax1 (*c2max, curh);
      if (cond1)
      {
	params->adx[i] = (FCELL) slp;
	params->ady[i] = (FCELL) oor;
	if (cond2)
	{
	  params->adxx[i] = (FCELL) curn;
	  params->adyy[i] = (FCELL) curh;
	  params->adxy[i] = (FCELL) curm;
	}
      }
    }				/* bmask == 1 */
  }
  return 1;
}