secpar2d.c

Go to the documentation of this file.

/*-
 * 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;
}