mymath.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include <InterViews/geometry.h>
#include "mymath.h"
#include "classreg.h"
#include "oc2iv.h"
#include <math.h>
#include <stdio.h>

extern int hoc_return_type_code;

static double distance_to_line(void*) {
    return MyMath::distance_to_line(*getarg(1), *getarg(2),
        *getarg(3), *getarg(4), *getarg(5), *getarg(6));
}
    
static double distance_to_line_segment(void*) {
    return MyMath::distance_to_line_segment(*getarg(1), *getarg(2),
        *getarg(3), *getarg(4), *getarg(5), *getarg(6));
}

static double inside(void*) {
    hoc_return_type_code = 2; // boolean
    return MyMath::inside(*getarg(1), *getarg(2), *getarg(3), *getarg(4),
        *getarg(5), *getarg(6));
}

int nrn_feround(int);
// return last rounding mode and set to given mode if 1,2,3,4.
// order is FE_DOWNWARD, FE_TONEAREST, FE_TOWARDZERO, FE_UPWARD
#if defined(HAVE_FENV_H) && defined(HAVE_FESETROUND)
#include <fenv.h>
static int round_mode[] = {FE_DOWNWARD, FE_TONEAREST, FE_TOWARDZERO, FE_UPWARD};
#endif
int nrn_feround(int mode) {
#if defined(HAVE_FENV_H) && defined(HAVE_FESETROUND)
    int oldmode = fegetround();
    int m;
    if (oldmode == FE_TONEAREST) {
        oldmode = 2;
    }else if (oldmode == FE_TOWARDZERO) {
        oldmode = 3;
    }else if (oldmode == FE_UPWARD) {
        oldmode = 4;
    }else if (oldmode == FE_DOWNWARD) {
        oldmode = 1;
    }else{
        assert(0);
    }
    if (mode > 0 && mode < 5) {
        nrn_assert(fesetround(round_mode[mode-1]) == 0);
    }
    return oldmode;
#else
    return 0;
#endif
}

static double feround(void*) {
    int arg = 0;
    hoc_return_type_code = 1; // integer
    if (ifarg(1)) {
        arg = (int)chkarg(1, 0, 4);
    }
    return (double) nrn_feround(arg);
}

static Member_func members[] = {
    "d2line", distance_to_line,
    "d2line_seg", distance_to_line_segment,
    "inside", inside,
    "feround", feround,
    0,0
};

static void* cons(Object*) {
    return NULL;
}

static void destruct(void*) {
}

void GUIMath_reg() {
    class2oc("GUIMath", cons, destruct, members, NULL, NULL, NULL);
}

double MyMath::anint(double x) {
    if (x < 0) {
        return ceil(x - .5);
    }else{
        return floor(x + .5);
    }
}
float MyMath::min(int count, const float* x) {
    float m=x[0];
    for (int i=1; i < count; ++i) {
        if (m > x[i]) {
            m = x[i];
        }
    }
    return m;
}

float MyMath::max(int count, const float* x) {
    float m=x[0];
    for (int i=1; i < count; ++i) {
        if (m < x[i]) {
            m = x[i];
        }
    }
    return m;
}
    
// within epsilon distance from the infinite line

bool MyMath::near_line(Coord x, Coord y,
   Coord x1, Coord y1, Coord x2, Coord y2,
   float epsilon)
{
//printf("near_line %g %g %g %g %g %g %g\n", x,y,x1,y1,x2,y2,epsilon);
    if (Math::equal(x, x1, epsilon) && Math::equal(y, y1, epsilon)) {
        return true;
    }
    if (Math::equal(x1, x2, epsilon) && Math::equal(y1, y2, epsilon)) {
        return false;
    }
    Coord d, norm, norm2, dot;
    Coord dx, dy, dx2, dy2;
    dx = x - x1;
    dy = y - y1;
    dx2 = x2 - x1;
    dy2 = y2 - y1;
//printf("%g %g %g %g\n", dx, dy, dx2, dy2);
    norm2 = dx2*dx2 + dy2*dy2;
    norm = dx*dx + dy*dy;
    dot = dx*dx2 + dy*dy2;
    d = norm - dot*dot/norm2;
//printf("near_line %g\n", d);
    if (d <= epsilon*epsilon) {
        return true;
    }else{
        return false;
    }
}

float MyMath::distance_to_line(Coord x, Coord y,
   Coord x1, Coord y1, Coord x2, Coord y2
  )
{
//printf("near_line %g %g %g %g %g %g %g\n", x,y,x1,y1,x2,y2,epsilon);
    Coord d, norm, norm2, dot;
    Coord dx, dy, dx2, dy2;
    dx = x - x1;
    dy = y - y1;
    dx2 = x2 - x1;
    dy2 = y2 - y1;
//printf("%g %g %g %g\n", dx, dy, dx2, dy2);
    norm2 = dx2*dx2 + dy2*dy2;
    norm = dx*dx + dy*dy;
    dot = dx*dx2 + dy*dy2;
    if (norm2 == 0) {norm2 = 1.;}
    d = norm - dot*dot/norm2;
    if (d < 0.) {
        return 0.;
    }
    return sqrt(d);
}

float MyMath::distance_to_line_segment(Coord x, Coord y,
   Coord x1, Coord y1, Coord x2, Coord y2
  )
{
    Coord norm, norm2, dot;
    Coord dx, dy, dx2, dy2;
    dx = x - x1;
    dy = y - y1;
    dx2 = x2 - x1;
    dy2 = y2 - y1;
    norm2 = dx2*dx2 + dy2*dy2;
    norm = dx*dx + dy*dy;

    if (norm2 == 0) {
        return sqrt(norm);
    }
    dot = dx*dx2 + dy*dy2;
    if (dot < 0) {
        return sqrt(norm);
    }else if (dot > norm2) {
        dx = x - x2;
        dy = y - y2;
        return sqrt(dx*dx + dy*dy);
    }else{
        dx = norm - dot*dot/norm2;
        if (dx <= 0) {
            return 0.;
        }
        return sqrt(dx);
    }
}

bool MyMath::near_line_segment(Coord x, Coord y,
   Coord x1, Coord y1, Coord x2, Coord y2,
   float epsilon)
{
    Coord l = x1, b = y1, r = x2, t = y2;
    MyMath::minmax(l, r);
    MyMath::minmax(b, t);
//printf("near_line_seg inside %d\n",inside(x, y, l-epsilon, b-epsilon, r+epsilon, t+epsilon));
//printf("near_line_seg near_line %d\n", near_line(x, y, l, b, r, t, epsilon));
    return 
        MyMath::inside(x, y, l-epsilon, b-epsilon, r+epsilon, t+epsilon)
        && MyMath::near_line(x, y, x1, y1, x2, y2, epsilon)
    ;
}
    
#if 0
void MyMath::round_range(Coord x1, Coord x2, double& y1, double& y2, int& ntic) {
    double d = x2 - x1;
    int e = 0;
    while (d > 10) {
        d /= 10.;
        ++e;
    }
    while (d < 10) {
        d *= 10.;
        --e;
    }

    if (d > 76.) {
        d = 100.; ntic = 5;
    }else if (d > 51.) {
        d = 75.; ntic = 3;
    }else if (d > 31.) {
        d = 50.; ntic = 5;
    }else if (d > 26.) {
        d = 30.; ntic = 3;
    }else if (d > 21.) {
        d = 25.; ntic = 5;
    }else if (d > 16) {
        d = 20; ntic = 4;
#if 0
    }else if (d > 11) {
        d = 15; ntic = 3;
    }else{
        d = 10.; ntic = 5;
    }
#else
    }else{
        int i = int(d);
        d = double(i);
        if (i > 6) {
            if ((i%2) == 0) {
                ntic = i/2;
            }else if ((i%3) == 0) {
                ntic = i/3;
            } else {
                i = i+1;
                d = double(i);
                ntic = i/2;
            }
        }else{
            ntic = i;
        }
    }
#endif
    
    while (e > 0) {
        d *= 10.;
        --e;
    }
    while (e < 0) {
        d /= 10.;
        ++e;
    }

    d /= double(ntic);
    y1 = d*double(int(x1/d - .5));
    y2 = d*double(int(x2/d + .5));
    ntic = int((y2 - y1)/d + .5);
}
#else
void MyMath::round_range(Coord x1, Coord x2, double& y1, double& y2, int& ntic) {
    double d = x2 - x1;
    d = pow(10, floor(log10(d)))/10;
    y1 = d*MyMath::anint(x1/d);
    y2 = d*MyMath::anint(x2/d);
    int i =  int((y2 - y1)/d + .5); // 10 < i < 100
//printf("%d %g\n", i, dx);
   for(;;) {
    if (i%3 == 0) {
        ntic = 3;
        return;
    }else if (i%4 == 0) {
        ntic = 4;
        return;
    }else if (i%5 == 0) {
        ntic = 5;
        return;
    }
    y1 -= d;
    y2 += d;
    i += 2;
   }
}
#endif

void MyMath::round_range_down(Coord x1, Coord x2, double& y1, double& y2, int& ntic) {
    double d = x2 - x1;
    double e = pow(10, floor(log10(d)))/10;
    int i =  int(d/e + .5);
    if (i > 20) {
        y1 = 5*e*ceil(x1/e/5 - .01);
        y2 = 5*e*floor(x2/e/5 + .01);
    }else{
        y1 = e*ceil(x1/e - .01);
        y2 = e*floor(x2/e + .01);
    }
    i = int((y2 - y1)/e + .5);
//printf("%d %g\n", i, e);
   for(;;) {
    if (i%3 == 0) {
        ntic = 3;
        return;
    }else if (i%4 == 0) {
        ntic = 4;
        return;
    }else if (i%5 == 0) {
        ntic = 5;
        return;
    }
    y1 -= e;
    ++i;
   }
}

double MyMath::round(float& x1, float& x2, int direction, int digits) {
    double d;
    if (x2 > x1) {
        d = x2 - x1;
    }else{
        d = Math::abs(x1);
    }
    double e = pow(10, floor(log10(d)) + 1 - digits);
    switch (direction) {
    case Expand:
        x1 = e*floor(x1/e);
        x2 = e*ceil(x2/e);
        break;
    case Contract:
        x1 = e*ceil(x1/e);
        x2 = e*floor(x2/e);
        break;
    case Lower:
        x1 = e*floor(x1/e);
        x2 = e*floor(x2/e);
        break;
    case Higher:
        x1 = e*ceil(x1/e);
        x2 = e*ceil(x2/e);
        break;
    }
    return e;
}

void MyMath::box(Requisition& req, Coord& x1, Coord& y1, Coord& x2, Coord& y2) {
    Requirement& rx = req.x_requirement();
    x1 = -rx.alignment() * rx.natural();
    x2 = x1 + rx.natural();
    Requirement& ry = req.y_requirement();
    y1 = -ry.alignment() * ry.natural();
    y2 = y1 + ry.natural();
}

bool MyMath::unit_normal(Coord x, Coord y, Coord* perp) {
    float d;
    d = sqrt(x*x + y*y);
    if (d < 1e-6) {
        perp[0] = 0.;
        perp[1] = 1.;
        return false;
    }
    perp[0] = y/d;
    perp[1] = -x/d;
    return true;
}