D7Dist.h

//
//  CD7Dist.c
//
//
//  Created by Herbert J. Bernstein and Lawrence C. Andrews on 3/04/17.
//
//

#ifndef D7DIST_H
#define D7DIST_H

#ifdef __cplusplus

extern "C" {

#endif




/* #define D7DIST_DEBUG */
#define D7DIST_NO_OUTER_PASS 

#include <math.h>
#include <float.h>

/*#define D7DIST_DEBUG */
#ifdef D7DIST_DEBUG
static int changed=0;
static double oldvalue;
#include <stdio.h>
#define report_double(prolog,value,epilog) \
oldvalue=value; fprintf(stderr,"%s%g%s",prolog,value,epilog);
#define report_integer(prolog,value,epilog) \
fprintf(stderr,"%s%d%s",prolog,value,epilog);
#define report_double_if_changed(prolog,value,epilog) \
changed=0; if (fabs(value-oldvalue)>1.e-8*(fabs(value)+fabs(oldvalue)+1.e-12)) {oldvalue=value; changed=1; fprintf(stderr,"%s%g%s",prolog,value,epilog);}
#define also_if_changed_report(prolog,value,epilog) \
if(changed) {fprintf(stderr,"%s%s%s",prolog,value,epilog);}
#define also_if_changed_report_integer(prolog,value,epilog) \
if(changed) {fprintf(stderr,"%s%d%s",prolog,value,epilog);}
#define also_if_changed_report_double(prolog,value,epilog) \
if(changed) {fprintf(stderr,"%s%g%s",prolog,value,epilog);}
#define also_if_changed_report_double_vector(prolog,value,epilog) \
if(changed) {fprintf(stderr,"%s[%g, %g, %g, %g, %g, %g, %g]%s",prolog,value[0],value[1],value[2],value[3],value[4],value[5],value[6],epilog);}
#define report_double_vector(prolog,value,epilog) \
{fprintf(stderr,"%s[%g, %g, %g, %g, %g, %g, %g]%s",prolog,value[0],value[1],value[2],value[3],value[4],value[5],value[6],epilog);}
#else
#define report_double(prolog,value,epilog)
#define report_integer(prolog,value,epilog)
#define report_double_if_changed(prolog,value,epilog)
#define also_if_changed_report(prolog,value,epilog)
#define also_if_changed_report_integer(prolog,value,epilog)
#define also_if_changed_report_double(prolog,value,epilog)
#define also_if_changed_report_double_vector(prolog,value,epilog)
#define report_double_vector(prolog,value,epilog)
#endif


#define CD7M_min(a,b) (((a)<(b))?(a):(b))
#define CD7M_max(a,b) (((a)<(b))?(b):(a))

#define ND7BND 9

#define D7P_1 0
#define D7P_2 1
#define D7P_3 2
#define D7P_4 3
#define D7P_5 4
#define D7P_6 5
#define D7P_7 6
#define D7P_8 7
#define D7P_9 8

#define ND7BDPRJ 9


/* D7Dist Boundary projectors */

static double d7prj[ND7BDPRJ][49]= {
    /*prj[D7P_1]  ,1,*/
    {0.3571428571428572,0.3571428571428572,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        0.3571428571428572,0.3571428571428572,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,0.8571428571428571,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        0.8571428571428571,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,0.8571428571428571,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,0.8571428571428571,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        0.8571428571428571},
    /*prj[D7P_2]   2,*/
    {0.8571428571428571,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,0.3571428571428572,0.3571428571428572,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,0.3571428571428572,0.357142571428572,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        0.8571428571428571,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,0.8571428571428571,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,0.8571428571428571,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        0.8571428571428571},
    /*prj[D7P_3]   3,*/
    {0.8571428571428571,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,0.8571428571428571,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,0.3571428571428572,
        0.3571428571428572,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,0.3571428571428572,
        0.3571428571428572,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,0.8571428571428571,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,0.8571428571428571,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        0.8571428571428571},
    /*prj[D7P_4]   4,*/
    {0.75,0.0,0.0,-0.25,0.0,0.25,0.25,
        0.0,0.6666666666666666,-0.3333333333333333,0.0,
        0.3333333333333333,0.0,0.0,
        0.0,-0.3333333333333333,0.6666666666666666,0.0,
        0.3333333333333333,0.0,0.0,-0.25,0.0,0.0,0.75,0.0,0.25,0.25,
        0.0,0.3333333333333333,0.3333333333333333,0.0,
        0.6666666666666666,0.0,0.0,0.25,0.0,0.0,0.25,0.0,0.75,-0.25,
        0.25,0.0,0.0,0.25,0.0,-0.25,0.75},
    /*prj[D7P_5]   5 */
    {0.6666666666666666,0.0,0.0,-0.3333333333333333,
        0.3333333333333333,0.0,0.0,0.0,0.75,-0.25,0.0,0.0,0.25,0.25,
        0.0,-0.25,0.75,0.0,0.0,0.25,0.25,
        -0.3333333333333333,0.0,0.0,0.6666666666666666,
        0.3333333333333333,0.0,0.0,
        0.3333333333333333,0.0,0.0,0.3333333333333333,
        0.6666666666666666,0.0,0.0,0.0,0.25,0.25,0.0,0.0,0.75,-0.25,
        0.0,0.25,0.25,0.0,0.0,-0.25,0.75},
    /*prj[D7P_6]   6 */
    {0.6666666666666666,0.0,-0.3333333333333333,0.0,0.0,
        0.3333333333333333,0.0,0.0,0.75,0.0,-0.25,0.25,0.0,0.25,
        -0.3333333333333333,0.0,0.6666666666666666,0.0,0.0,
        0.3333333333333333,0.0,0.0,-0.25,0.0,0.75,0.25,0.0,0.25,
        0.0,0.25,0.0,0.25,0.75,0.0,-0.25,
        0.3333333333333333,0.0,0.3333333333333333,0.0,0.0,
        0.6666666666666666,0.0,0.0,0.25,0.0,0.25,-0.25,0.0,0.75},
    /*prj[D7P_7]   7 */
    {0.75,0.0,-0.25,0.0,0.25,0.0,0.25,
        0.0,0.6666666666666666,0.0,-0.3333333333333333,0.0,
        0.3333333333333333,0.0,-0.25,0.0,0.75,0.0,0.25,0.0,0.25,
        0.0,-0.3333333333333333,0.0,0.6666666666666666,0.0,
        0.3333333333333333,0.0,0.25,0.0,0.25,0.0,0.75,0.0,-0.25,
        0.0,0.3333333333333333,0.0,0.3333333333333333,0.0,
        0.6666666666666666,0.0,0.25,0.0,0.25,0.0,-0.25,0.0,0.75},
    /*prj[D7P_8]   8 */
    {0.6666666666666666,-0.3333333333333333,0.0,0.0,0.0,0.0,
        0.3333333333333333,
        -0.3333333333333333,0.6666666666666666,0.0,0.0,0.0,0.0,
        0.3333333333333333,0.0,0.0,0.75,-0.25,0.25,0.25,0.0,
        0.0,0.0,-0.25,0.75,0.25,0.25,0.0,
        0.0,0.0,0.25,0.25,0.75,-0.25,0.0,
        0.0,0.0,0.25,0.25,-0.25,0.75,0.0,
        0.3333333333333333,0.3333333333333333,0.0,0.0,0.0,0.0,
        0.6666666666666666},
    /*prj[D7P_9]   9 */
    {0.75,-0.25,0.0,0.0,0.25,0.25,0.0,
        -0.25,0.75,0.0,0.0,0.25,0.25,0.0,
        0.0,0.0,0.6666666666666666,-0.3333333333333333,0.0,0.0,
        0.3333333333333333,
        0.0,0.0,-0.3333333333333333,0.6666666666666666,0.0,0.0,
        0.3333333333333333,0.25,0.25,0.0,0.0,0.75,-0.25,0.0,
        0.25,0.25,0.0,0.0,-0.25,0.75,0.0,
        0.0,0.0,0.3333333333333333,0.3333333333333333,0.0,0.0,
        0.6666666666666666}};

static double d7prj_perp[ND7BDPRJ][49]= {
    /*prj_perp[D7P_1]   1 */
    {0.6428571428571429,-0.3571428571428572,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        -0.3571428571428572,0.6428571428571429,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428},
    /*prj_perp[D7P_2]   2 */
    {0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.6428571428571429,-0.3571428571428572,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,-0.3571428571428572,0.6428571428571429,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428},
    /*prj_perp[D7P_3]   3 */
    {0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,0.6428571428571429,
        -0.3571428571428572,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        0.1428571428571428,0.1428571428571428,-0.3571428571428572,
        0.6428571428571429,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428,
        -0.1428571428571428,-0.1428571428571428,-0.1428571428571428,
        -0.1428571428571428,0.1428571428571428,0.1428571428571428,
        0.1428571428571428},
    /*prj_perp[D7P_4]   4 */
    {0.25,0.0,0.0,0.25,0.0,-0.25,-0.25,
        0.0,0.3333333333333333,0.3333333333333333,0.0,
        -0.3333333333333333,0.0,0.0,
        0.0,0.3333333333333333,0.3333333333333333,0.0,
        -0.3333333333333333,0.0,0.0,
        0.25,0.0,0.0,0.25,0.0,-0.25,-0.25,
        0.0,-0.3333333333333333,-0.3333333333333333,0.0,
        0.3333333333333333,0.0,0.0,
        -0.25,0.0,0.0,-0.25,0.0,0.25,0.25,
        -0.25,0.0,0.0,-0.25,0.0,0.25,0.25},
    /*prj_perp[D7P_5]   5 */
    {0.3333333333333333,0.0,0.0,0.3333333333333333,
        -0.3333333333333333,0.0,0.0,
        0.0,0.25,0.25,0.0,0.0,-0.25,-0.25,
        0.0,0.25,0.25,0.0,0.0,-0.25,-0.25,
        0.3333333333333333,0.0,0.0,0.3333333333333333,
        -0.3333333333333333,0.0,0.0,
        -0.3333333333333333,0.0,0.0,-0.3333333333333333,
        0.3333333333333333,0.0,0.0,
        0.0,-0.25,-0.25,0.0,0.0,0.25,0.25,
        0.0,-0.25,-0.25,0.0,0.0,0.25,0.25},
    /*prj_perp[D7P_6]   6 */
    {0.3333333333333333,0.0,0.3333333333333333,0.0,0.0,
        -0.3333333333333333,0.0,0.0,0.25,0.0,0.25,-0.25,0.0,-0.25,
        0.3333333333333333,0.0,0.3333333333333333,0.0,0.0,
        -0.3333333333333333,0.0,0.0,0.25,0.0,0.25,-0.25,0.0,-0.25,
        0.0,-0.25,0.0,-0.25,0.25,0.0,0.25,
        -0.3333333333333333,0.0,-0.3333333333333333,0.0,0.0,
        0.3333333333333333,0.0,0.0,-0.25,0.0,-0.25,0.25,0.0,0.25},
    /*prj_perp[D7P_7]   7 */
    {0.25,0.0,0.25,0.0,-0.25,0.0,-0.25,
        0.0,0.3333333333333333,0.0,0.3333333333333333,0.0,
        -0.3333333333333333,0.0,0.25,0.0,0.25,0.0,-0.25,0.0,-0.25,
        0.0,0.3333333333333333,0.0,0.3333333333333333,0.0,
        -0.3333333333333333,0.0,-0.25,0.0,-0.25,0.0,0.25,0.0,0.25,
        0.0,-0.3333333333333333,0.0,-0.3333333333333333,0.0,
        0.3333333333333333,0.0,-0.25,0.0,-0.25,0.0,0.25,0.0,0.25},
    /*prj_perp[D7P_8]   8 */
    {0.3333333333333333,0.3333333333333333,0.0,0.0,0.0,0.0,
        -0.3333333333333333,
        0.3333333333333333,0.3333333333333333,0.0,0.0,0.0,0.0,
        -0.3333333333333333,0.0,0.0,0.25,0.25,-0.25,-0.25,0.0,
        0.0,0.0,0.25,0.25,-0.25,-0.25,0.0,
        0.0,0.0,-0.25,-0.25,0.25,0.25,0.0,
        0.0,0.0,-0.25,-0.25,0.25,0.25,0.0,
        -0.3333333333333333,-0.3333333333333333,0.0,0.0,0.0,0.0,
        0.3333333333333333},
    /*prj_perp[D7P_9]   9 */
    {0.25,0.25,0.0,0.0,-0.25,-0.25,0.0,
        0.25,0.25,0.0,0.0,-0.25,-0.25,0.0,
        0.0,0.0,0.3333333333333333,0.3333333333333333,0.0,0.0,
        -0.3333333333333333,
        0.0,0.0,0.3333333333333333,0.3333333333333333,0.0,0.0,
        -0.3333333333333333,-0.25,-0.25,0.0,0.0,0.25,0.25,0.0,
        -0.25,-0.25,0.0,0.0,0.25,0.25,0.0,
        0.0,0.0,-0.3333333333333333,-0.3333333333333333,0.0,0.0,
        0.3333333333333333}};



/* The following matrices are the transformation
 matrices that may be applied at the associated
 boundaries.  *** An addition post permutation may
 be needed to preserve D7 ordering *** */

/* D7 post permuations {1,2,3,4,5,6,7},   Ident
                       {1,2,4,3,6,5,7},   D7M_3
                       {1,3,2,4,5,7,6},   D7M_2
                       {1,3,4,2,7,5,6},   D7M_1.D7M_2
                       {1,4,2,3,6,7,5},
                       {1,4,3,2,7,6,5},
                       {3,1,2,4,7,5,6},
                       {4,1,2,3,7,6,5}
 */

#define D7Refl_1 0
#define D7Refl_2 6
#define D7Refl_3 2
#define D7Refl_4 1
#define D7Refl_term -1

static int D7Rord[24]= {0,1,2,3,4,5,12,18,6,7,8,9,10,11,13,14,15,16,17,19,20,21,22,23};


static int D7Refl[24][49]={
    /*D7Refl_1 {1,2,3,4,5,6,7},   Ident*/
    /*DR_0*/ {
        1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1},
    /*D7Refl_4 {1,2,4,3,6,5,7},   c <-> d */
    /*DR_1*/ {
        1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1},
    /*D7Refl_3 {1,3,2,4,5,7,6},   b <-> c */
    /*DR_2*/ {
        1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0},
    /* {1,3,4,2,7,5,6} */
    /*DR_3*/ {
        1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0},
    /* {1,4,2,3,6,7,5} */
    /*DR_4*/ {
        1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0},
    /* {1,4,3,2,6.7,5} */
    /*DR_5*/ {
        1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,1,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0},
    /*D7Refl_2 {2,1,3,4,6,5,7},   a <-> b */
    /*DR_6*/ {
        0,1,0,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1},
    /*DR_7*/ {
        0,1,0,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1},
    /*DR_8*/ {
        0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0},
    /*DR_9*/ {
        0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0},
    /*DR_10*/ {
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0},
    /*DR_11*/ {
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,1,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0},
    /* {3,1,2,4,7,5,6} */
    /*DR_12*/ {
        0,0,1,0,0,0,0,
        1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0},
    /*DR_13*/ {
        0,0,1,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0},
    /*DR_14*/ {
        0,0,1,0,0,0,0,
        0,1,0,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0},
    /*DR_15*/ {
        0,0,1,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0},
    /*DR_16*/ {
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1}, 
    /*DR_17*/ {
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        0,1,0,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1}, 
    /*DR_18*/ {
        0,0,0,1,0,0,0,
        1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0}, 
    /*DR_19*/ {
        0,0,0,1,0,0,0,
        1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0}, 
    /*DR_20*/ {
        0,0,0,1,0,0,0,
        0,1,0,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0}, 
    /*DR_21*/ {
        0,0,0,1,0,0,0,
        0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0}, 
    /*DR_22*/ {
        0,0,0,1,0,0,0,
        0,0,1,0,0,0,0,
        1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1}, 
    /*DR_23*/ {
        0,0,0,1,0,0,0,
        0,0,1,0,0,0,0,
        0,1,0,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1}
};


/*  The 24 D7 permuations, the first 8 of which are seen in the
    boundary transforms */

static int D7Perm[24][8] = {
    
    /* D7Perm_1 {1,2,3,4,5,6,7},   Ident   OK */
    {D7Refl_1,D7Refl_term,0,0,0,0,0,0},
    /* D7Perm_2 {1,2,4,3,6,5,7}            OK */
    {D7Refl_4,D7Refl_term,0,0,0,0,0,0},
    /* D7Perm_3 {1,3,2,4,5,7,6}            OK */
    {D7Refl_3,D7Refl_term,0,0,0,0,0,0},
    /* D7Perm_4 {1,3,4,2,7,5,6},           OK */
    {D7Refl_4,D7Refl_3,D7Refl_term,0,0,0,0,0},
    /* D7Perm_5 {1,4,2,3,6,7,5}            OK */
    {D7Refl_3,D7Refl_4,D7Refl_term,0,0,0,0,0},
    /* D7Perm_6  {1,4,3,2,7,6,5}           OK */
    {D7Refl_4,D7Refl_3,D7Refl_4,D7Refl_term,0,0,0,0},
    /* D7Perm_7  {3,1,2,4,7,5,6}           OK */
    {D7Refl_2,D7Refl_3,D7Refl_term,0,0,0,0,0},
    /* D7Perm_8  {4,1,2,3,7,6,5}           OK */
    {D7Refl_2,D7Refl_3,D7Refl_4,D7Refl_term,0,0,0,0},
    /* D7Perm_9  {2, 1, 3, 4}              OK */
    {D7Refl_2,D7Refl_term,0,0,0,0,0,0},
    /* D7Perm_10 {2, 1 ,4, 3}              OK */
    {D7Refl_4,D7Refl_2,D7Refl_term,0,0,0,0,0},
    /* D7Perm_11 {2, 3, 1, 4}              OK */
    {D7Refl_3,D7Refl_2,D7Refl_term,0,0,0,0,0},
    /* D7Perm_12 {2, 3, 4, 1}              OK */
    {D7Refl_4,D7Refl_3,D7Refl_2,D7Refl_term,0,0,0,0},
    /* D7Perm_13 {2, 4, 1, 3}              OK */
    {D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_term,0,0,0,0},
    /* D7Perm_14 {2, 4, 3, 1}              OK */
    {D7Refl_4,D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_term,0,0,0},
    /* D7Perm_15 {3, 1, 4, 2}              OK */
    {D7Refl_4,D7Refl_2,D7Refl_3,D7Refl_term,0,0,0,0},
    /* D7Perm_16 {3, 2, 1, 4}              OK */
    {D7Refl_3,D7Refl_2,D7Refl_3,D7Refl_term,0,0,0,0},
    /* D7Perm_17 {3, 2, 4, 1}              OK */
    {D7Refl_4,D7Refl_3,D7Refl_2,D7Refl_3,D7Refl_term,0,0,0},
    /* D7Perm_18 {3, 4, 1, 2}              OK */
    {D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_3,D7Refl_term,0,0,0},
    /* D7Perm_19 {3, 4, 2, 1}              OK */
    {D7Refl_4,D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_3,D7Refl_term,0,0},
    /* D7Perm_20 {4, 1, 3, 2}              OK */
    {D7Refl_2,D7Refl_4,D7Refl_3,D7Refl_4,D7Refl_term,0,0,0},
    /* D7Perm_21 {4, 2, 1, 3}              OK */
    {D7Refl_2,D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_term,0,0,0},
    /* D7Perm_22 {4, 2, 3, 1}              OK */
    {D7Refl_2,D7Refl_4,D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_term,0,0},
    /* D7Perm_23 {4, 3, 1, 2}              OK */
    {D7Refl_2,D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_3,D7Refl_term,0,0},
    /* D7Perm_24 {4, 3, 2, 1}              OK */
    {D7Refl_4,D7Refl_2,D7Refl_3,D7Refl_4,D7Refl_2,D7Refl_3,D7Refl_term,0}
};

static int D7MS[9][49] = {
    
    /* D7M_1 d1 <-> d2, d5 <-> d6 {2,1,3,4,6,5,7} */
    {0,1,0,0,0,0,0,
        1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1},

    /* D7M_2 d2 <-> d3, d6 <-> d7 {1,3,2,4,5,7,6} */
    {1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1,
        0,0,0,0,0,1,0},

    /* D7M_3 d3 <-> d4, d5 <-> d6 {1,2,4,3,6,5,7}*/
    {1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1},
    
    /* D7M_4 d1->d2, d2->d3, d3->d6, d4->d7, d5->d1, d6->d4, d7->2d2+2d3-d5
       and post permutations 1,2,3,4,5,6,7,8
     */
    {0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1,
        1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,2,2,0,-1,0,0},
    
    /* D7M_5 d2->d4, d3->d5, d4->d1, d5->d2, d6->d3, d7->-2d2-2d3+d5+2d6+2d7
       and post permutations 1,2,3,4,5,6,7,8 */
    {1,0,0,0,0,0,0,
        0,0,0,1,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,0,0,1,
        0,1,0,0,0,0,0,
        0,0,1,0,0,0,0,
        2,0,0,2,-1,0,0},
    
    /* D7M_6 d2->d3, d3->d5, d4->d7, d5->d2, d6->d4, d6->2d1+2d3-d6
       and post permutations 1,2,3,4,5,6,7,8 */
    
    {1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        2,0,2,0,0,-1,0},

    /* D7M_7 d2->d3, d3->d5, d4->d1, d5->d2, d6->d4, d6-2d1+2d3-d6
     and post permutations 1,2,3,4,5,6,7,8, all are case 6 cells */
    
    {1,0,0,0,0,0,0,
        0,0,1,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,0,1,
        0,1,0,0,0,0,0,
        0,0,0,1,0,0,0,
        2,0,2,0,0,-1,0},
    
    /* D7M_8 d3->d5, d4->d6, d5->d3, d6->d4, d7->2d1+2d2-d7
     and post permutations 1,2,3,5,8,7 */
    
    {1,0,0,0,0,0,0,
        0,1,0,0,0,0,0,
        0,0,0,0,1,0,0,
        0,0,0,0,0,1,0,
        0,0,1,0,0,0,0,
        0,0,0,1,0,0,0,
        2,2,0,0,0,0,-1},
    
    /* D7M_9 d1->d3, d2->d6, d3->d5, d5->-2d1-2d2+2d5+2d6+d7,d6->d2,d7->d1 */
    
    {0,0,1,0,0,0,0,
        0,0,0,0,0,1,0,
        0,0,0,0,1,0,0,
        0,0,0,1,0,0,0,
        0,0,2,2,0,0,-1,
        0,1,0,0,0,0,0,
        1,0,0,0,0,0,0}

};




#undef M_E
#define M_1 0
#define M_2 1
#define M_3 2
#define M_4 3
#define M_5 4
#define M_6 5
#define M_7 6
#define M_8 7
#define M_9 8
#define M_A 9
#define M_B 10
#define M_C 11
#define M_D 12
#define M_E 13
#define M_F 14


#define NMPGS            65
#define NPGS             18



/* Test if outside D7 region
   return 0 if outside */

int d7test(double g[7]);


/* Compute the dot product of 2 D7 vectors */

static double d7dotprod(double v1[7], double v2[7]) {
    
    return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]+v1[3]*v2[3]+v1[4]*v2[4]+v1[5]*v2[5]+v1[6]*v2[6];
    
}


/* Compute the simple Euclidean distance between 2 D7 vectors */

static double d7eucldistsq(double v1[7], double v2[7]) {
    
    double vtemp;
    int ii;
    double distsq;
    
    distsq = 0.;
    
    for (ii = 0; ii < 7; ii++ ) {
        vtemp = v1[ii]-v2[ii];
        distsq += vtemp*vtemp;
    }
    return distsq;
}

static double d7eucldist(double v1[7], double v2[6]) {
    
    return sqrt(d7eucldistsq(v1,v2));
    
}



#define CD7M_min2(a,b) (((a)<(b))?(a):(b))
#define CD7M_min3(a,b,c) ((CD7M_min2(a,b)<(c))?(CD7M_min2(a,b)):(c))
#define CD7M_min4(a,b,c,d) ((CD7M_min3(a,b,c)<(d))?(CD7M_min3(a,b,c)):(d))
#define CD7M_min5(a,b,c,d,e) ((CD7M_min4(a,b,c,d)<(e))?(CD7M_min4(a,b,c,d)):(e))
#define CD7M_min6(a,b,c,d,e,f) ((CD7M_min5(a,b,c,d,e)<(f))?(CD7M_min5(a,b,c,d,e)):(f))
#define CD7M_min7(a,b,c,d,e,f,g) ((CD7M_min6(a,b,c,d,e,f)<(g))?(CD7M_min6(a,b,c,d,e,f)):(g))
#define CD7M_min8(a,b,c,d,e,f,g,h) ((CD7M_min7(a,b,c,d,e,f,g)<(h))?(CD7M_min7(a,b,c,d,e,f,g)):(h))


/* Macro versions of distances */

#define CD7M_normsq(v) v[0]*v[0]+v[1]*v[1]+v[2]*v[2]+v[3]*v[3]+v[4]*v[4]+v[5]*v[5]
#define CD7M_norm(v) sqrt(CD7M_normsq(v))
#define CD7M_d7normsq(v) v[0]*v[0]+v[1]*v[1]+v[2]*v[2]+v[3]*v[3]+v[4]*v[4]+v[5]*v[5]+v[6]*v[6]
#define CD7M_d7norm(v) sqrt(CD7M_d7normsq(v))

#define CD7M_gdistsq(v1,v2) \
( \
(v1[0]-v2[0])*(v1[0]-v2[0])+\
(v1[1]-v2[1])*(v1[1]-v2[1])+\
(v1[2]-v2[2])*(v1[2]-v2[2])+\
(v1[3]-v2[3])*(v1[3]-v2[3])+\
(v1[4]-v2[4])*(v1[4]-v2[4])+\
(v1[5]-v2[5])*(v1[5]-v2[5]))

#define CD7M_gdist(v1,v2) sqrt(CD7M_gdistsq(v1,v2))

#define CD7M_d7eucldistsq(v1,v2) \
( \
(v1[0]-v2[0])*(v1[0]-v2[0])+\
(v1[1]-v2[1])*(v1[1]-v2[1])+\
(v1[2]-v2[2])*(v1[2]-v2[2])+\
(v1[3]-v2[3])*(v1[3]-v2[3])+\
(v1[4]-v2[4])*(v1[4]-v2[4])+\
(v1[5]-v2[5])*(v1[5]-v2[5])+\
(v1[6]-v2[6])*(v1[6]-v2[6])\
)

#define CD7M_d7eucldist(v1,v2) sqrt(CD7M_d7eucldistsq(v1,v2))






/*     Compute the best distance between 2 d6 vectors
 allowing for permulations of d1, d2, d3, d4
 */

#define CD7M_d7eucldistsq_byelem(v11,v12,v13,v14,v15,v16,v17,v21,v22,v23,v24,v25,v26,v27) \
fabs((v11-v21)*(v11-v21)+(v12-v22)*(v12-v22)+(v13-v23)*(v13-v23) + \
(v14-v24)*(v14-v24)+(v15-v25)*(v15-v25)+(v16-v26)*(v16-v26) + (v17-v27)*(v17-v27))

#define CD7M_d7prods_byelem(v11,v12,v13,v14,v15,v16,v17,v21,v22,v23,v24,v25,v26,v27) \
2.*(v11*v21+v12*v22+v13*v23+v14*v24+v15*v25+v16*v26+v17*v27)

#define CD7M_d71234distsq(v1,v2) \
CD7M_min3(\
CD7M_min8(\
CD7M_d7eucldistsq_byelem(v1[0],v1[1],v1[2],v1[3],v1[4],v1[5],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[0],v1[1],v1[3],v1[2],v1[5],v1[4],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[0],v1[2],v1[1],v1[3],v1[4],v1[6],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[0],v1[2],v1[3],v1[1],v1[6],v1[4],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[0],v1[3],v1[1],v1[2],v1[5],v1[6],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[0],v1[3],v1[2],v1[1],v1[6],v1[5],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[2],v1[0],v1[1],v1[3],v1[6],v1[4],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[3],v1[0],v1[1],v1[2],v1[6],v1[5],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6])),\
CD7M_min8(\
CD7M_d7eucldistsq_byelem(v1[1],v1[0],v1[2],v1[3],v1[5],v1[4],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[1],v1[0],v1[3],v1[2],v1[4],v1[5],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[1],v1[2],v1[0],v1[3],v1[5],v1[6],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[1],v1[2],v1[3],v1[0],v1[6],v1[5],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[1],v1[3],v1[0],v1[2],v1[4],v1[6],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[1],v1[3],v1[2],v1[0],v1[6],v1[4],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[2],v1[0],v1[3],v1[1],v1[4],v1[6],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[2],v1[1],v1[0],v1[3],v1[6],v1[5],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6])),\
CD7M_min8(\
CD7M_d7eucldistsq_byelem(v1[2],v1[1],v1[3],v1[0],v1[5],v1[6],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[2],v1[3],v1[0],v1[1],v1[4],v1[5],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[2],v1[3],v1[1],v1[0],v1[5],v1[4],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[3],v1[0],v1[2],v1[1],v1[5],v1[6],v1[4],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[3],v1[1],v1[0],v1[2],v1[6],v1[4],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[3],v1[1],v1[2],v1[0],v1[4],v1[6],v1[5],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[3],v1[2],v1[0],v1[1],v1[5],v1[4],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]),\
CD7M_d7eucldistsq_byelem(v1[3],v1[2],v1[1],v1[0],v1[4],v1[5],v1[6],\
v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6])))

#define CD7M_d71234dist(v1,v2) sqrt(CD7M_d71234distsq(v1,v2))

static double d71234distsq(double v1[7], double v2[7]) {
    double distsq;
    double dwnby;
    double dtrial[24];
    int i;
    distsq = v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2] + v1[3]*v1[3] + v1[4]*v1[4] + v1[5]*v1[5] + v1[6]*v1[6]
    + v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2] + v2[3]*v2[3] + v2[4]*v2[4] + v2[5]*v2[5] + v2[6]*v2[6];
    dtrial[0] = CD7M_d7prods_byelem(v1[0],v1[1],v1[2],v1[3],v1[4],v1[5],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[1] = CD7M_d7prods_byelem(v1[0],v1[1],v1[3],v1[2],v1[5],v1[4],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[2] = CD7M_d7prods_byelem(v1[0],v1[2],v1[1],v1[3],v1[4],v1[6],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[3] = CD7M_d7prods_byelem(v1[0],v1[2],v1[3],v1[1],v1[6],v1[4],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[4] = CD7M_d7prods_byelem(v1[0],v1[3],v1[1],v1[2],v1[5],v1[6],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[5] = CD7M_d7prods_byelem(v1[0],v1[3],v1[2],v1[1],v1[6],v1[5],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[6] = CD7M_d7prods_byelem(v1[1],v1[0],v1[2],v1[3],v1[5],v1[4],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[7] = CD7M_d7prods_byelem(v1[1],v1[0],v1[3],v1[2],v1[4],v1[5],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[8] = CD7M_d7prods_byelem(v1[1],v1[2],v1[0],v1[3],v1[5],v1[6],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[9] = CD7M_d7prods_byelem(v1[1],v1[2],v1[3],v1[0],v1[6],v1[5],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[10] = CD7M_d7prods_byelem(v1[1],v1[3],v1[0],v1[2],v1[4],v1[6],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[11] = CD7M_d7prods_byelem(v1[1],v1[3],v1[2],v1[0],v1[6],v1[4],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[12] = CD7M_d7prods_byelem(v1[2],v1[0],v1[1],v1[3],v1[6],v1[4],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[13] = CD7M_d7prods_byelem(v1[2],v1[0],v1[3],v1[1],v1[4],v1[6],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[14] = CD7M_d7prods_byelem(v1[2],v1[1],v1[0],v1[3],v1[6],v1[5],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[15] = CD7M_d7prods_byelem(v1[2],v1[1],v1[3],v1[0],v1[5],v1[6],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[16] = CD7M_d7prods_byelem(v1[2],v1[3],v1[0],v1[1],v1[4],v1[5],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[17] = CD7M_d7prods_byelem(v1[2],v1[3],v1[1],v1[0],v1[5],v1[4],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[18] = CD7M_d7prods_byelem(v1[3],v1[0],v1[1],v1[2],v1[6],v1[5],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[19] = CD7M_d7prods_byelem(v1[3],v1[0],v1[2],v1[1],v1[5],v1[6],v1[4],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[20] = CD7M_d7prods_byelem(v1[3],v1[1],v1[0],v1[2],v1[6],v1[4],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[21] = CD7M_d7prods_byelem(v1[3],v1[1],v1[2],v1[0],v1[4],v1[6],v1[5],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[22] = CD7M_d7prods_byelem(v1[3],v1[2],v1[0],v1[1],v1[5],v1[4],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    dtrial[23] = CD7M_d7prods_byelem(v1[3],v1[2],v1[1],v1[0],v1[4],v1[5],v1[6],v2[0],v2[1],v2[2],v2[3],v2[4],v2[5],v2[6]);
    
    dwnby = dtrial[0];
    for (i = 1; i < 24; i++) if (dtrial[i] > dwnby) dwnby = dtrial[i];
    return (fabs(distsq - dwnby));
}

#define d71234dist(v1,v2) sqrt(d71234distsq(v1,v2))


/*     Compute the best distance between 2 G6 vectors
 allowing for permulations of g1, g2, g3 as
 well as sign changes
 */


static void d7cpyvn(int n, double src[], double dst[] ) {
    int i;
    for (i = 0; i < n; i++) {
        dst[i] = src[i];
    }
}


static void imv7 (double v1[7], int m[49], double v2[7]) {
    int i, j;
    double sum;
    for (i = 0; i < 7; i++) {
        sum = 0.0;
        for(j=0; j < 7; j++) {
            sum = sum + ((double)m[7*i+j])*v1[j];
        }
        v2[i] = sum;
    }
}


static void rmv7 (double v1[7], double m[49], double v2[7]) {
    int i, j;
    double sum;
    for (i = 0; i < 7; i++) {
        sum = 0.0;
        for(j=0; j < 7; j++) {
            sum = sum + m[7*i+j]*v1[j];
        }
        v2[i] = sum;
    }
}




static void d7twoPminusI(double pg[7], double g[7], double gout[7]) {
    int i;
    for (i=0; i<7; i++){
        gout[i] = 2.*pg[i]-g[i];
    }
}


/*     Map a G6 vector onto the boundaries after
       applying the 24-way unfolding */

#ifdef D7DIST_NO_OUTER_PASS
#define NREFL_OUTER_FULL 8
#define NREFL_OUTER_MIN 1
#define NREFL_OUTER_MID 4
#else
#define NREFL_OUTER_FULL 24
#define NREFL_OUTER_MID 16
#define NREFL_OUTER_MIN 8
#endif

static double d7bddist(double gvec[7],int bdnum) {
    
    /* inward normals in D7 */
    double D7N_1x[7]={-1./sqrt(2.),1./sqrt(2.),0.,0.,0.,0.,0.};
    double D7N_2x[7]={0.,-1./sqrt(2.),1./sqrt(2.),0.,0.,0.,0.};
    double D7N_3x[7]={0.,0.,-1./sqrt(2.),1./sqrt(2.),0.,0.,0.};
    double D7N_4x[7]={-3/(2*sqrt(21)),2/sqrt(21),2/sqrt(21),-3/(2*sqrt(21)),
                       -2/sqrt(21),3/(2*sqrt(21)),3/(2*sqrt(21))};
    double D7N_5x[7]={2/sqrt(21),-3/(2*sqrt(21)),-3/(2*sqrt(21)),2/sqrt(21),
                       -2/sqrt(21),3/(2*sqrt(21)),3/(2*sqrt(21))};
    double D7N_6x[7]={2/sqrt(21),-3/(2*sqrt(21)),2/sqrt(21),-3/(2*sqrt(21)),
                       3/(2*sqrt(21)),-2/sqrt(21),3/(2*sqrt(21))};
    double D7N_7x[7]={-3/(2*sqrt(21)),2/sqrt(21),-3/(2*sqrt(21)),2/sqrt(21),
                       3/(2*sqrt(21)),-2/sqrt(21),3/(2*sqrt(21))};
    double D7N_8x[7]={2/sqrt(21),2/sqrt(21),-3/(2*sqrt(21)),-3/(2*sqrt(21)),
                       3/(2*sqrt(21)),3/(2*sqrt(21)),-2/sqrt(21)};
    double D7N_9x[7]={-3/(2*sqrt(21)),-3/(2*sqrt(21)),2/sqrt(21),2/sqrt(21),
                       3/(2*sqrt(21)),3/(2*sqrt(21)),-2/sqrt(21)};
    
    
    if (bdnum < ND7BND) {
        
        switch(bdnum) {
            case(D7P_1): return -d7dotprod(D7N_1x,gvec); break;
            case(D7P_2): return -d7dotprod(D7N_2x,gvec); break;
            case(D7P_3): return -d7dotprod(D7N_3x,gvec); break;
            case(D7P_4): return -d7dotprod(D7N_4x,gvec); break;
            case(D7P_5): return -d7dotprod(D7N_5x,gvec); break;
            case(D7P_6): return -d7dotprod(D7N_6x,gvec); break;
            case(D7P_7): return -d7dotprod(D7N_7x,gvec); break;
            case(D7P_8): return -d7dotprod(D7N_8x,gvec); break;
            case(D7P_9): return -d7dotprod(D7N_9x,gvec); break;
                
            default: return DBL_MAX; break;
        }
        
    }
    
    return DBL_MAX;
}




/* Get the minimum distance to the boundaries */

static double d7minbddist(double gvec[7]) {
    int ii;
    double dists[9];
    double minbd;
    dists[0] = fabs(d7bddist(gvec,D7P_1));
    dists[1] = fabs(d7bddist(gvec,D7P_2));
    dists[2] = fabs(d7bddist(gvec,D7P_3));
    dists[3] = fabs(d7bddist(gvec,D7P_4));
    dists[4] = fabs(d7bddist(gvec,D7P_5));
    dists[5] = fabs(d7bddist(gvec,D7P_6));
    dists[6] = fabs(d7bddist(gvec,D7P_7));
    dists[7] = fabs(d7bddist(gvec,D7P_8));
    dists[8] = fabs(d7bddist(gvec,D7P_9));
    
    minbd = dists[0];
    for (ii=1; ii<9; ii++) {
        minbd = CD7M_min(minbd,dists[ii]);
    }
    return minbd;
}



static void d7bdmaps(double gvec[7],
            double dists[ND7BND],
            int iord[ND7BND],
            double pgs[ND7BND][7],
            double rgs[ND7BND][7],
            double mpgs[ND7BND][7],
            double mvecs[ND7BND][7],
            double maxdist,
            int * ngood) {
    
    int jj, itemp, igap, idone;
    
    /* double D7N_1x[7]={-19./sqrt(1106),-5./sqrt(1106),-12./sqrt(1106),-12./sqrt(1106),12./sqrt(1106),12./sqrt(1106),12./sqrt(1106)};
    double D7N_2x[7]={-12./sqrt(1106),-19./sqrt(1106),-5./sqrt(1106),-12./sqrt(1106),12./sqrt(1106),12./sqrt(1106),12./sqrt(1106)};
    double D7N_3x[7]={-12./sqrt(1106),-12./sqrt(1106),-19./sqrt(1106),-5./sqrt(1106),12./sqrt(1106),12./sqrt(1106),12./sqrt(1106)};
    double D7N_4x[7]={-3./sqrt(39),1./sqrt(39),1./sqrt(39),-3./sqrt(39),-1./sqrt(39),3./sqrt(39),3./sqrt(39)};
    double D7N_5x[7]={1./sqrt(39),-3./sqrt(39),-3./sqrt(39),1./sqrt(39),-1./sqrt(39),3./sqrt(39),3./sqrt(39)};
    double D7N_6x[7]={0,-1./2.,0,-1./2.,1./2.,0,1./2};
    double D7N_7x[7]={-3./sqrt(39),1./sqrt(39),-3./sqrt(39),1./sqrt(39),3./sqrt(39),-1./sqrt(39),3./sqrt(39)};
    double D7N_8x[7]={-2./4.,-2./4.,-1./4.,-1./4.,1./4.,1./4.,2./4.};
    double D7N_9x[7]={-3./sqrt(39),-3./sqrt(39),1./sqrt(39),1./sqrt(39),3./sqrt(39),3./sqrt(39),-1./sqrt(39)}; */
    
    /* inward normals in D7 */
    double D7N_1x[7]={-1./sqrt(2.),1./sqrt(2.),0.,0.,0.,0.,0.};
    double D7N_2x[7]={0.,-1./sqrt(2.),1./sqrt(2.),0.,0.,0.,0.};
    double D7N_3x[7]={0.,0.,-1./sqrt(2.),1./sqrt(2.),0.,0.,0.};
    double D7N_4x[7]={-3/(2*sqrt(21)),2/sqrt(21),2/sqrt(21),-3/(2*sqrt(21)),
                       -2/sqrt(21),3/(2*sqrt(21)),3/(2*sqrt(21))};
    double D7N_5x[7]={2/sqrt(21),-3/(2*sqrt(21)),-3/(2*sqrt(21)),2/sqrt(21),
                       -2/sqrt(21),3/(2*sqrt(21)),3/(2*sqrt(21))};
    double D7N_6x[7]={2/sqrt(21),-3/(2*sqrt(21)),2/sqrt(21),-3/(2*sqrt(21)),
                       3/(2*sqrt(21)),-2/sqrt(21),3/(2*sqrt(21))};
    double D7N_7x[7]={-3/(2*sqrt(21)),2/sqrt(21),-3/(2*sqrt(21)),2/sqrt(21),
                       3/(2*sqrt(21)),-2/sqrt(21),3/(2*sqrt(21))};
    double D7N_8x[7]={2/sqrt(21),2/sqrt(21),-3/(2*sqrt(21)),-3/(2*sqrt(21)),
                       3/(2*sqrt(21)),3/(2*sqrt(21)),-2/sqrt(21)};
    double D7N_9x[7]={-3/(2*sqrt(21)),-3/(2*sqrt(21)),2/sqrt(21),2/sqrt(21),
                       3/(2*sqrt(21)),3/(2*sqrt(21)),-2/sqrt(21)};
    
    
    dists[D7P_1] = fabs(d7dotprod(D7N_1x,gvec));
    dists[D7P_2] = fabs(d7dotprod(D7N_2x,gvec));
    dists[D7P_3] = fabs(d7dotprod(D7N_3x,gvec));
    dists[D7P_4] = fabs(d7dotprod(D7N_4x,gvec));
    dists[D7P_5] = fabs(d7dotprod(D7N_5x,gvec));
    dists[D7P_6] = fabs(d7dotprod(D7N_6x,gvec));
    dists[D7P_7] = fabs(d7dotprod(D7N_7x,gvec));
    dists[D7P_8] = fabs(d7dotprod(D7N_8x,gvec));
    dists[D7P_9] = fabs(d7dotprod(D7N_9x,gvec));
                
    
    for (jj=0; jj < ND7BND; jj++) {
        iord[jj] = jj;
    }

    
    igap = ND7BND;
    while (igap > 1) {
        igap = igap/2;
        idone = 0;
        while (!idone) {
            idone = 1;
            for (jj=0; jj < ND7BND-igap; jj+=igap) {
                if (dists[iord[jj]] > dists[iord[jj+igap]]) {
                    idone = 0;
                    itemp = iord[jj];
                    iord[jj] = iord[jj+igap];
                    iord[jj+igap] = itemp;
                }
            }
        }
    }

    *ngood = ND7BND;
    for (jj = 0; jj < ND7BND; jj++ ) {
        rmv7(gvec, d7prj[jj], pgs[jj]);
        d7twoPminusI(pgs[jj],gvec,rgs[jj]);
        imv7(pgs[jj], D7MS[jj], mpgs[jj]);
        imv7(gvec, D7MS[jj],mvecs[jj]);
        if (dists[jj] > maxdist) (*ngood)--;
    }
    
}

/* revised D7Dist_2bds
 
   Compute lengths of paths from g_lft to g_rgt through bd_up and bd_dwn
 
 
 */

static double D7Dist_2bds_rev(double g_lft[7], double g_rgt[7],
                           double pg_lft_up[7], double mpg_lft_up[7],
                           double pg_rgt_up[7], double mpg_rgt_up[7], int bd_up,
                           double pg_rgt_dwn[7], double mpg_rgt_dwn[7],
                           double pg_lft_dwn[7], double mpg_lft_dwn[7], int bd_dwn,
                           double dist) {
    
    double bddist_lft_up, bddist_lft_dwn, bddist_right_up, bddist_right_dwn;
    double alpha_lft_up, alpha_lft_dwn;
    double bdint_lft_up[7], bdint_lft_dwn[7], mbdint_lft_up[7], mbdint_lft_dwn[7];
    int ii;
    double gvec_lft_bd_lft_up, gvec_lft_mbd_lft_up, gvec_lft_bd_lft_dwn, gvec_lft_mbd_lft_dwn;
    double gvec_lft_bd_rgt_dwn, gvec_lft_mbd_rgt_dwn;
    double gvec_rgt_bd_lft_up, gvec_rgt_mbd_lft_up, gvec_rgt_bd_lft_dwn, gvec_rgt_mbd_lft_dwn;
    double gvec_rgt_bd_rgt_dwn, gvec_rgt_mbd_rgt_dwn;
    double gvec_lft_bd_lft_upbd_rgt_dwn, gvec_lft_bd_lft_upmbd_rgt_dwn;
    double gvec_lft_mbd_lft_upbd_rgt_dwn, gvec_lft_mbd_lft_upmbd_rgt_dwn;
    double bd_lft_upbd_rgt_dwn, bd_lft_upmbd_rgt_dwn;
    double mbd_lft_upbd_rgt_dwn, mbd_lft_upmbd_rgt_dwn;

    double dist_glft_bd_up, dist_grgt_bd_up;
    double dist_glft_bd_dwn, dist_grgt_bd_dwn;
    double dist_bd_up_bd_dwn, dist_bd_up_mbd_dwn, dist_mbd_up_bd_dwn, dist_mbd_up_mbd_dwn;
    
    
    /* distances from g_lft and g_rgt to bd_up and bd_up */
    
    bddist_lft_up = fabs(d7bddist(g_lft, bd_up));
    bddist_lft_dwn = fabs(d7bddist(g_lft, bd_dwn));
    bddist_right_up = fabs(d7bddist(g_rgt, bd_up));
    bddist_right_dwn = fabs(d7bddist(g_rgt, bd_dwn));
    
    /* the portions of a minimal mirror path off each boundary */
    
    if (bddist_lft_up+bddist_right_up < bddist_lft_up*1.e-10) {
        alpha_lft_up = 0.;
    } else {
        alpha_lft_up = bddist_lft_up/(bddist_lft_up+bddist_right_up);
    }
    if (bddist_lft_dwn+bddist_right_dwn < bddist_lft_dwn*1.e-10) {
        alpha_lft_dwn = 0.;
    } else {
        alpha_lft_dwn = bddist_lft_dwn/(bddist_lft_dwn+bddist_right_dwn);
    }
    
    if (alpha_lft_up > 1.) alpha_lft_up = 1.;
    if (alpha_lft_dwn > 1.) alpha_lft_dwn = 1.;
    
    
    if (alpha_lft_up < 0.) alpha_lft_up = 0.;
    if (alpha_lft_dwn < 0.) alpha_lft_dwn = 0.;
    
    for (ii = 0; ii < 7; ii++) {
        bdint_lft_up[ii] = pg_lft_up[ii] + alpha_lft_up*(pg_rgt_up[ii] - pg_lft_up[ii]);
        bdint_lft_dwn[ii] = pg_lft_dwn[ii] + alpha_lft_dwn*(pg_rgt_dwn[ii] - pg_lft_dwn[ii]);
    }
    imv7(bdint_lft_up, D7MS[bd_up], mbdint_lft_up);
    imv7(bdint_lft_dwn, D7MS[bd_dwn], mbdint_lft_dwn);
    
    /* possible routes
     g_lft - bdint_lft_up - g_rgt
     g_lft - mbdint_lft_up - g_rgt
     g_lft - bdint_lft_dwn - g_rgt
     g_lft - mbdint_lft_dwn - g_rgt
     
     min(g_lft - mbdint_lft_up, bdint_lft_up)+min(bdint_left_up - mdint_right_dn - g_rgt
     g_lft - mbdint_lft_up - mbdint21 - g_rgt
     g_lft - mbdint_lft_up - bdint22 - g_rgt
     g_lft - mbdint_lft_up - mbdint22 - g_rgt
     
     */
     
     
    gvec_lft_bd_lft_up = d71234dist(g_lft,bdint_lft_up);
    gvec_lft_mbd_lft_up = d71234dist(g_lft,mbdint_lft_up);
    gvec_lft_bd_lft_dwn = d71234dist(g_lft,bdint_lft_dwn);
    gvec_lft_mbd_lft_dwn = d71234dist(g_lft,mbdint_lft_dwn);
    gvec_rgt_bd_lft_up = d71234dist(g_rgt,bdint_lft_up);
    gvec_rgt_mbd_lft_up = d71234dist(g_rgt,mbdint_lft_up);
    gvec_rgt_bd_lft_dwn = d71234dist(g_rgt,bdint_lft_dwn);
    gvec_rgt_mbd_lft_dwn = d71234dist(g_rgt,mbdint_lft_dwn);
    
    dist_bd_up_bd_dwn = d71234dist(bdint_lft_up,bdint_lft_dwn);
    dist_bd_up_mbd_dwn = d71234dist(bdint_lft_up,mbdint_lft_dwn);
    if (dist_bd_up_mbd_dwn < dist_bd_up_bd_dwn ) dist_bd_up_bd_dwn = dist_bd_up_mbd_dwn;
    dist_mbd_up_bd_dwn = d71234dist(mbdint_lft_up,bdint_lft_dwn);
    if (dist_mbd_up_bd_dwn < dist_bd_up_bd_dwn ) dist_bd_up_bd_dwn = dist_mbd_up_bd_dwn;
    dist_mbd_up_mbd_dwn = d71234dist(mbdint_lft_up,mbdint_lft_dwn);
    if (dist_mbd_up_mbd_dwn < dist_bd_up_bd_dwn ) dist_bd_up_bd_dwn = dist_mbd_up_mbd_dwn;
    
    
    if (gvec_lft_bd_lft_up + gvec_rgt_bd_lft_up < dist) dist = gvec_lft_bd_lft_up + gvec_rgt_bd_lft_up;
    if (gvec_lft_mbd_lft_up + gvec_rgt_mbd_lft_up < dist) dist = gvec_lft_mbd_lft_up + gvec_rgt_mbd_lft_up;
    if (gvec_lft_mbd_lft_up + gvec_rgt_bd_lft_up < dist) dist = gvec_lft_mbd_lft_up + gvec_rgt_bd_lft_up;
    if (gvec_lft_mbd_lft_up + gvec_rgt_mbd_lft_up < dist) dist = gvec_lft_mbd_lft_up + gvec_rgt_mbd_lft_up;
    
    if (gvec_lft_bd_lft_dwn + gvec_rgt_bd_lft_dwn < dist) dist = gvec_lft_bd_lft_dwn + gvec_rgt_bd_lft_dwn;
    if (gvec_lft_mbd_lft_dwn + gvec_rgt_mbd_lft_dwn < dist) dist = gvec_lft_mbd_lft_dwn + gvec_rgt_mbd_lft_dwn;
    if (gvec_lft_mbd_lft_dwn + gvec_rgt_bd_lft_dwn < dist) dist = gvec_lft_mbd_lft_dwn + gvec_rgt_bd_lft_dwn;
    if (gvec_lft_mbd_lft_dwn + gvec_rgt_mbd_lft_dwn < dist) dist = gvec_lft_mbd_lft_dwn + gvec_rgt_mbd_lft_dwn;
    
    dist_glft_bd_up = gvec_lft_bd_lft_up;
    if (gvec_lft_mbd_lft_up < dist_glft_bd_up) dist_glft_bd_up = gvec_lft_mbd_lft_up;
    dist_glft_bd_dwn = gvec_lft_bd_lft_dwn;
    if (gvec_lft_mbd_lft_dwn < dist_glft_bd_dwn) dist_glft_bd_dwn = gvec_lft_mbd_lft_dwn;
    dist_grgt_bd_up = gvec_rgt_bd_lft_up;
    if (gvec_rgt_mbd_lft_up < dist_grgt_bd_up) dist_grgt_bd_up = gvec_rgt_mbd_lft_up;
    dist_grgt_bd_dwn = gvec_rgt_bd_lft_dwn;
    if (gvec_rgt_mbd_lft_dwn < dist_grgt_bd_dwn) dist_grgt_bd_dwn = gvec_rgt_mbd_lft_dwn;
    
    if (gvec_lft_mbd_lft_up + dist_bd_up_bd_dwn + dist_grgt_bd_dwn < dist ) dist = gvec_lft_mbd_lft_up + dist_bd_up_bd_dwn + dist_grgt_bd_dwn;
    if (gvec_lft_mbd_lft_dwn + dist_bd_up_bd_dwn + dist_grgt_bd_up < dist ) dist = gvec_lft_mbd_lft_dwn + dist_bd_up_bd_dwn + dist_grgt_bd_up;
    
    
    
    return dist;
    
}




#include <stdio.h>

/* Compute the minimal distance between gvec1 and gvec2 going
 through bd1 and bd2, assuming rgvec1 and rgevc2 are their
 reflections in those boundaries.
 
 */



#define DCUT 0.999995
#define fudge(d) DCUT*d


/*
 Compute tminimal distance between two Delaunay-reduced
 vectors in the Delaunay Cone following the embedding paths
 to the 9 boundaries
 */



static double D7Dist_pass(double gvec1[7],double gvec2[7],double dist) {
    double dists1[ND7BND];
    double pgs1[ND7BND][7], rgs1[ND7BND][7], mpgs1[ND7BND][7], mvecs1[ND7BND][7];
    double dists2[ND7BND];
    double pgs2[ND7BND][7], rgs2[ND7BND][7], mpgs2[ND7BND][7], mvecs2[ND7BND][7];
    int iord1[ND7BND],iord2[ND7BND];
    double mindists1;
    double mindists2;
    int jx1,jx2;
    int j1,j2;
    int ngood1,ngood2;
    double maxdist;
    
    maxdist = fudge(dist);
    
    d7bdmaps(gvec1,dists1,iord1,pgs1,rgs1,mpgs1,mvecs1,maxdist,&ngood1);
    d7bdmaps(gvec2,dists2,iord2,pgs2,rgs2,mpgs2,mvecs2,maxdist,&ngood2);
    
    mindists1 = d7minbddist(gvec1);
    mindists2 = d7minbddist(gvec2);
    
    if (mindists1+mindists2 > dist) return dist;
    
    if (mindists1+mindists2 < maxdist) {
        for (jx1 = 0; jx1 < ngood1; jx1++) {
            double d1;
            j1 = iord1[jx1];
            if (j1 < D7P_4 || j1 ==D7P_7 || j1 ==D7P_9) continue;
            d1 = dists1[j1];

            if (d1 < maxdist){
                dist = CD7M_min(dist,d71234dist(gvec2,mpgs1[j1])+d1);
            }
        }
        for (jx2 = 0; jx2 < ngood2; jx2++) {
            double d2;
            j2 = iord2[jx2];
            if (j2 < D7P_4 || j2 ==D7P_7 || j2 ==D7P_9) continue;
            d2 = dists2[j2];
            if (d2 < maxdist) {
                dist = CD7M_min(dist,(d71234dist(gvec1,mpgs2[j2])+d2));
                
            }
        }
    }
    
    maxdist = fudge(dist);
    for (jx1 = 0; jx1 < ngood1; jx1++) {
        double d1;
        j1 = iord1[jx1];
        if (j1 < D7P_4 || j1 ==D7P_7 || j1 ==D7P_9 ) continue;
        d1 = dists1[j1];
        if (d1 < maxdist) {
            for (jx2 = 0; jx2 < ngood2; jx2++) {
                double d2;
                j2 = iord2[jx2];
                if (j2 < D7P_4 || j2 ==D7P_7 || j2 ==D7P_9 ) continue;
                d2 = dists2[j2];
                if (d2 < maxdist) {
                    dist = D7Dist_2bds_rev(gvec1, gvec2, pgs1[j1], mpgs1[j1], pgs2[j1], mpgs2[j1],j1,
                                        pgs2[j2], mpgs2[j2], pgs1[j2], mpgs1[j2], j2, dist);
                    
                }
            }
        }
    }
    return dist;
}


/*
     Compute the minimal distance between two Delaunay-reduced
     vectors in the D7 Cone following the embedding paths
     to the 7 major boundaries
 */





double D7Dist(double * gvec1,double * gvec2);


#ifdef __cplusplus

}

#endif
#endif /*D7DIST_H */