 /*
  * Khoros: $Id: lvslope.c,v 1.1 1991/05/10 15:43:20 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvslope.c,v 1.1 1991/05/10 15:43:20 khoros Exp $";
#endif

 /*
  * $Log: lvslope.c,v $
 * Revision 1.1  1991/05/10  15:43:20  khoros
 * Initial revision
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvslope.c
 >>>>
 >>>>      Program Name: vslope
 >>>>
 >>>> Date Last Updated: Thu Mar  7 15:34:12 1991 
 >>>>
 >>>>          Routines: lvslope - the library call for vslope
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define radians_to_degrees 57.295780
#define two_pi 6.2831853
#define pi 3.1415926
#define pi_over_2 1.5707963
#define pi_over_12 0.26179939
#define rad_eq_8_deg 0.1396263
#define max_float 1.0e+20
#define min_float (1/max_float)

#define CHECKLIBINPUT(program, img) {                                         \
        if (!proper_map_enable (program, img, VFF_MAP_OPTIONAL, FALSE)) {     \
            fprintf (stderr, "error in lvslope: input image may not have a ");\
            fprintf (stderr, "forced map\n");                                 \
            return (0);                                                       \
        }                                                                     \
        if (!proper_num_images (program, img, 1, FALSE)) {                    \
            fprintf (stderr, "error in lvslope: input image must contain ");  \
            fprintf (stderr, "only 1 image\n");                               \
            return (0);                                                       \
        }                                                                     \
        if (!proper_loc_type (program, img, VFF_LOC_IMPLICIT, FALSE)) {       \
            fprintf (stderr, "error in lvslope: input image may not contain");\
            fprintf (stderr, " explicit location data\n");                    \
            return (0);                                                       \
        }                                                                     \
        if ((img->data_storage_type != VFF_TYP_1_BYTE) &&                     \
            (img->data_storage_type != VFF_TYP_2_BYTE) &&                     \
            (img->data_storage_type != VFF_TYP_4_BYTE) &&                     \
            (img->data_storage_type != VFF_TYP_FLOAT)) {                      \
            fprintf (stderr, "error in lvslope: input data type must be ");   \
            fprintf (stderr, "byte, short, int, or float\n");                 \
            return (0);                                                       \
        }                                                                     \
        if (img->pixsizx <= 0) {                                              \
            fprintf (stderr,                                                  \
                "error in lvslope: x pixel size must be greater than 0\n");   \
            return (0);                                                       \
        }                                                                     \
        if (img->pixsizy <= 0) {                                              \
            fprintf (stderr,                                                  \
                "error in lvslope: y pixel size must be greater than 0\n");   \
            return (0);                                                       \
        }                                                                     \
        if (img->pixsizx != img->pixsizy) {                                   \
            fprintf (stderr,                                                  \
                "error in lvslope: x and y pixel sizes must be the same\n");  \
            return (0);                                                       \
        }                                                                     \
}

static void compute_normal();
static void find_norm_ul();
static void find_norm_ur();
static void find_norm_lr();
static void find_norm_ll();
static void crossp();
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvslope - library call for vslope
*
* Purpose:
*    
*    Compute slope and aspect images from elevation data.
*    
*    
* Input:
*    
*         1. img1  -  A pointer to the input viff structure.
*    
*         2. slope_flag  -  When set to 1, this flag causes the  slope
*         image to be created.  0 means no slope image is created.
*    
*         3. slope_option  -  Determines the form of the slope output:
*         0=degrees, 1=radians, 2=percent.
*    
*         4. aspect_flag  -  When set  to  1,  this  flag  causes  the
*         aspect  image  to  be  created.   0 means no aspect image is
*         created.
*    
*         5. aspect_option  -  Determines the form of the aspect  out-
*         put: 0=degrees, 1=radians, 2=quadrants.
*    
*         6. no_aspect  -  Value given for the aspect of a  flat  sur-
*         face.   The  aspect gives the direction the slope is facing.
*         When there is no slope the, aspect is given  this  value  to
*         denote that it is undefined.
*    
*    
* Output:
*    
*         1. img2  -  A double pointer to the  slope  viff  structure.
*         Only used if slope_flag is set.
*    
*         2. img3  -  A double pointer to the aspect  viff  structure.
*         Only used if aspect flag is set.
*    
*    
*
* Written By: Per Lysne
*    
*    
****************************************************************/


/* -library_def */
int
lvslope(img1, slope_flag, slope_option, aspect_flag, aspect_option, no_aspect,
        img2, img3)
struct xvimage *img1;
int slope_flag, slope_option, aspect_flag, aspect_option;
float no_aspect;
struct xvimage **img2, **img3;
/* -library_def_end */

/* -library_code */
{
    /*
     * GEOMETRY ASSUMPTIONS:
     *
     * 1. X axis increases from left to right.
     * 2. Y axis increases from back to front.
     * 3. Z axis increases from bottom to top.
     *
     * Result: left handed coordinate system.
     */

    char *program = "lvslope";
    int i, j, tmp_int;
    int band, num_data_bands;
    int rows, cols;
    float a, b, c, rise, run;
    float *surface, *slope, *aspect, pixsiz;
    double angle;
    struct xvimage *slope_img, *aspect_img;

    /*
     * Check the inputs.
     */
    if (img1 == NULL) {
        fprintf (stderr, "error in lvslope: bad input image\n");
        return (0);
    }

    if ((slope_flag != 0) && (slope_flag != 1)) {
        fprintf (stderr, "error in lvslope: bad slope flag\n");
        return (0);
    }

    if ((slope_option != 0) && (slope_option != 1) && (slope_option != 2)) {
        fprintf (stderr, "error in lvslope: bad slope option\n");
        return (0);
    }

    if ((aspect_flag != 0) && (aspect_flag != 1)) {
        fprintf (stderr, "error in lvslope: bad aspect flag\n");
        return (0);
    }

    if ((aspect_option != 0) && (aspect_option != 1) && (aspect_option != 2)) {
        fprintf (stderr, "error in lvslope: bad aspect option\n");
        return (0);
    }

    CHECKLIBINPUT (program, img1);

    /*
     * Initialize
     */
    num_data_bands = img1->num_data_bands;
    cols = img1->row_size;
    rows = img1->col_size;
    pixsiz = img1->pixsizx;
    surface = (float *) (img1->imagedata);

    /*
     * Allocate and zero the space for the slope data.
     */
    slope = (float *) malloc (num_data_bands*rows*cols*sizeof(float));
    if (slope == NULL) {
        fprintf (stderr, "error in lvslope: ");
        fprintf (stderr, "malloc failed for slope data\n");
        return (0);
    }
    bzero (slope,num_data_bands*rows*cols*sizeof(float));

    /*
     * Allocate and zero the space for the aspect data.
     */
    aspect = (float *) malloc (num_data_bands*rows*cols*sizeof(float));
    if (aspect == NULL) {
        fprintf (stderr, "error in lvslope: ");
        fprintf (stderr, "malloc failed for aspect image\n");
        return (0);
    }
    bzero (aspect,num_data_bands*rows*cols*sizeof(float));

    /*
     * Loop for each band in the input image.
     */
    for (band=0; band<num_data_bands; band++) {
        /*
         * Loop through each pixel in the elevation image and find the
         * corresponding slope and aspect for that point.
         */
        for (i=0; i<rows; i++) {
            for (j=0; j<cols; j++) {

                compute_normal(surface,j,i,rows,cols,band,pixsiz,&a,&b,&c);

                /*
                 * Calculate the slope. 
                 */
                run = sqrt (a*a + b*b);
                rise = c;

                /* Slope is in degrees. */
                if (slope_option == 0) { 
                    if (run < min_float) {
                        slope[BPIXEL(band,j,i,rows,cols)] = 0.0;
                    }
                    else {
                        angle =  atan (rise/run);
                        slope[BPIXEL(band,j,i,rows,cols)] = 
                              90.0 - (radians_to_degrees * (float) angle);
                    }
                }

                /* Slope is in radians. */
                else if (slope_option == 1) {
                    if (run < min_float) {
                        slope[BPIXEL(band,j,i,rows,cols)] = 0.0;
                    }
                    else {
                        angle =  atan (rise/run);
                        slope[BPIXEL(band,j,i,rows,cols)] =
                              pi_over_2 - (float) angle;
                    }
                }

                /* Slope is in percent rise. */
                else {
                    if (rise < min_float) {
                        slope[BPIXEL(band,j,i,rows,cols)] = max_float;
                    }
                    else {
                        /*
                         * Remember: precent rise of surface, not normal vector,
                         * so run & rise are reversed! 
                         */
                        slope[BPIXEL(band,j,i,rows,cols)] = (run/rise) * 100.0;
                    }
                }

                /*
                 * Calculate the aspect.
                 */
                /* Aspect is in degrees. */
                if (aspect_option == 0) {

                    /* No aspect, surface is flat. */
                    if ((fabs(a) < min_float) && (fabs(b) < min_float)) {
                        aspect[BPIXEL(band,j,i,rows,cols)] = (float) no_aspect;
                    }

                    /* Aspect angle is 90 or 270 degress - atan is undefined. */
                    else if (fabs(a) < min_float) {
                        if (b >= 0.0) {
                            aspect[BPIXEL(band,j,i,rows,cols)] = 90.0;
                        }
                        else {
                            aspect[BPIXEL(band,j,i,rows,cols)] = 270.0;
                        }
                    }

                    /* Aspect can be found using atan. */
                    else {
                        /* Tan is defined from -90 to 90. */
                        angle = atan (b/a);
                        if (a>=0 && b>=0) {  /* 1st quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] =
                                            radians_to_degrees * (float) angle;
                        }
                        else if (a<0 && b>=0) {  /* 2nd quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] =
                                  180.0 + (radians_to_degrees * (float) angle);
                        }
                        else if (a<0 && b<0) {  /* 3rd quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] =
                                  180.0 + (radians_to_degrees * (float) angle);
                        }
                        else { /* 4th quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] =
                                  360.0 + (radians_to_degrees * (float) angle);
                        }
                    }
                }

                /* Aspect is in radians. */
                else if (aspect_option == 1) {

                    /* No aspect, surface is flat. */
                    if ((fabs(a) < min_float) && (fabs(b) < min_float)) {
                        aspect[BPIXEL(band,j,i,rows,cols)] = (float) no_aspect;
                    }
 
                    /*Aspect angle is pi/2 or 3pi/2 degress-atan is undefined.*/
                    else if (fabs(a) < min_float) {
                        if (b >= 0.0) {
                            aspect[BPIXEL(band,j,i,rows,cols)] = pi_over_2;
                        }
                        else {
                            aspect[BPIXEL(band,j,i,rows,cols)] = 3.0*pi_over_2;
                        }
                    }

                    /* Aspect can be found using atan. */
                    else {
                        /* Tan is defined from -pi/2 to pi/2. */
                        angle = atan (b/a);
                        if (a>=0 && b>=0) {  /* 1st quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] = (float) angle;
                        }
                        else if (a<0 && b>=0) {  /* 2nd quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] = 
                                   pi + (float)angle;
                        }
                        else if (a<0 && b<0) {  /* 3rd quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] =
                                   pi + (float)angle;
                        }
                        else { /* 4th quadrant */
                            aspect[BPIXEL(band,j,i,rows,cols)] =
                                   two_pi + (float)angle;
                        }
                    }
                }

                /* Aspect is in quadrants. */
                else {

                    /* No aspect, surface is flat. */
                    if ((fabs(a) < min_float) && (fabs(b) < min_float)) {
                        aspect[BPIXEL(band,j,i,rows,cols)] = (float) no_aspect;
                    }

                    /* Aspect angle is 90 or 270 degress - atan is undefined. */
                    else if (fabs(a) < min_float) {
                        if (b >= 0.0) {
                            aspect[BPIXEL(band,j,i,rows,cols)] = 7.0;
                        }
                        else {
                            aspect[BPIXEL(band,j,i,rows,cols)] = 19.0;
                        }
                    }

                    /* Aspect can be found using atan. */
                    else {

                        /* Tan is defined from -pi/2 to pi/2. */
                        angle = atan (b/a);
                        /* In the 1st quadrant, angle=angle. */
                        if (a<0 && b>=0) {  /* 2nd quadrant */
                            angle = pi + angle;
                        }
                        else if (a<0 && b<0) {  /* 3rd quadrant */
                            angle = pi + angle;
                        }
                        else if (a>=0 && b<0) { /* 4th quadrant */
                            angle = two_pi + angle;
                        }

                        /*
                         * Map the angle into the proper quadrant. Remember that
                         * I am misusing the term and there are actually 24
                         * quadrants.
                         */
                        angle = angle + rad_eq_8_deg;
                        if (angle <= two_pi) {
                            tmp_int = (int) (angle/pi_over_12);
                        }
                        else {
                            tmp_int = 0;
                        }
                        aspect[BPIXEL(band,j,i,rows,cols)] = 
                               (float) tmp_int + 1.0;
                    }
                }
            }
        }
    } /* end for loop */

    if (slope_flag) {
        /*
         * Create a new viff images to return the slope data.
         */
        slope_img = createimage (rows, cols, VFF_TYP_FLOAT, 1,
                                 num_data_bands,"slope image created by vslope",
                                 img1->map_row_size, img1->map_col_size,
                                 img1->map_scheme, img1->map_storage_type, 
                                 VFF_LOC_IMPLICIT, 0);
        if (slope_img->maps != NULL) free (slope_img->maps);
        if (slope_img->imagedata != NULL) free (slope_img->imagedata);
        slope_img->maps = (char *) img1->maps;
        slope_img->imagedata = (char *) slope;
        *img2 = slope_img;
    }

    if (aspect_flag) {
        /*
         * Create a new viff images to return the aspect data.
         */
        aspect_img = createimage (rows, cols, VFF_TYP_FLOAT, 1,
                                num_data_bands,"aspect image created by vslope",
                                  img1->map_row_size, img1->map_col_size,
                                  img1->map_scheme, img1->map_storage_type,
                                  VFF_LOC_IMPLICIT, 0);
        if (aspect_img->maps != NULL) free (aspect_img->maps);
        if (aspect_img->imagedata != NULL) free (aspect_img->imagedata);
        aspect_img->maps = (char *) img1->maps;
        aspect_img->imagedata = (char *) aspect;
        *img3 = aspect_img;
    }

    return(1);
}



static void compute_normal(surface,x,y,rows,cols,band,pixsiz,a,b,c)
float *surface,pixsiz,*a,*b,*c;
int x,y,rows,cols;
{
    /*
       Compute the surface normal at each point in SURFACE and store
       this vector in the vector image NORMAL at the corresponding
       location. 
    */
    float norm1[3],norm2[3],norm3[3],norm4[3];

    if (x == 0 && y == 0)
      {
        /* Upper left corner */
        find_norm_ul(surface[BPIXEL(band,x,y+1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x+1,y,rows,cols)],pixsiz,norm1);
        *a = norm1[0];
        *b = norm1[1];
        *c = norm1[2];
      }
    else if (x == cols-1 && y == 0)
      {
        /* Upper right corner */
        find_norm_ur(surface[BPIXEL(band,x-1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y+1,rows,cols)],pixsiz,norm1);
        *a = norm1[0];
        *b = norm1[1];
        *c = norm1[2];
      }
    else if (x == cols-1 && y == rows-1)
      {
        /* Lower right corner */
        find_norm_lr(surface[BPIXEL(band,x,y-1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x-1,y,rows,cols)],pixsiz,norm1);
        *a = norm1[0];
        *b = norm1[1];
        *c = norm1[2];
      }
    else if (x == 0 && y == rows-1)
      {
        /* Lower left corner */
        find_norm_ll(surface[BPIXEL(band,x+1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y-1,rows,cols)],pixsiz,norm1);
        *a = norm1[0];
        *b = norm1[1];
        *c = norm1[2];
      }
    else if (y == 0)
      {
        /* Top edge */
        find_norm_ul(surface[BPIXEL(band,x,y+1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x+1,y,rows,cols)],pixsiz,norm1);
        find_norm_ur(surface[BPIXEL(band,x-1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y+1,rows,cols)],pixsiz,norm2);
        *a = 0.5*(norm1[0]+norm2[0]);
        *b = 0.5*(norm1[1]+norm2[1]);
        *c = 0.5*(norm1[2]+norm2[2]);
      }
    else if (y == rows-1)
      {
        /* Bottom edge */
        find_norm_lr(surface[BPIXEL(band,x,y-1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x-1,y,rows,cols)],pixsiz,norm1);
        find_norm_ll(surface[BPIXEL(band,x+1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y-1,rows,cols)],pixsiz,norm2);
        *a = 0.5*(norm1[0]+norm2[0]);
        *b = 0.5*(norm1[1]+norm2[1]);
        *c = 0.5*(norm1[2]+norm2[2]);
      }
    else if (x == 0)
      {
        /* Left edge */
        find_norm_ul(surface[BPIXEL(band,x,y+1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x+1,y,rows,cols)],pixsiz,norm1);
        find_norm_ll(surface[BPIXEL(band,x+1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y-1,rows,cols)],pixsiz,norm2);
        *a = 0.5*(norm1[0]+norm2[0]);
        *b = 0.5*(norm1[1]+norm2[1]);
        *c = 0.5*(norm1[2]+norm2[2]);
      }
    else if (x == cols-1)
      {
        /* Right edge */
        find_norm_ur(surface[BPIXEL(band,x-1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y+1,rows,cols)],pixsiz,norm1);
        find_norm_lr(surface[BPIXEL(band,x,y-1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x-1,y,rows,cols)],pixsiz,norm2);
        *a = 0.5*(norm1[0]+norm2[0]);
        *b = 0.5*(norm1[1]+norm2[1]);
        *c = 0.5*(norm1[2]+norm2[2]);
      }
    else
      {
        /* In the middle somewhere */
        find_norm_ul(surface[BPIXEL(band,x,y+1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x+1,y,rows,cols)],pixsiz,norm1);
        find_norm_ur(surface[BPIXEL(band,x-1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y+1,rows,cols)],pixsiz,norm2);
        find_norm_lr(surface[BPIXEL(band,x,y-1,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x-1,y,rows,cols)],pixsiz,norm3);
        find_norm_ll(surface[BPIXEL(band,x+1,y,rows,cols)],
                  surface[BPIXEL(band,x,y,rows,cols)],
                  surface[BPIXEL(band,x,y-1,rows,cols)],pixsiz,norm4);
        *a = 0.25*(norm1[0]+norm2[0]+norm3[0]+norm4[0]);
        *b = 0.25*(norm1[1]+norm2[1]+norm3[1]+norm4[1]);
        *c = 0.25*(norm1[2]+norm2[2]+norm3[2]+norm4[2]);
      }

    return;
}



static void find_norm_ul(a,b,c,size,norm)
float a,b,c,size,*norm;
{
    /* Compute normal for upper left corner geometry */
    float v1[3],v2[3];

    /* Load the vectors */
    v1[0] = 0;      /* X component */
    v1[1] = size;   /* Y component */
    v1[2] = a-b;    /* Z component */
    v2[0] = size;   /* X component */
    v2[1] = 0;      /* Y component */
    v2[2] = c-b;    /* Z component */

    /* Obtain cross-product */
    crossp(v1,v2,norm);

    return;
} 



static void find_norm_ur(a,b,c,size,norm)
float a,b,c,size,*norm;
{
    /* Compute normal for upper right corner geometry */
    float v1[3],v2[3];

    /* Load the vectors */
    v1[0] = -size;   /* X component */
    v1[1] = 0;      /* Y component */
    v1[2] = a-b;    /* Z component */
    v2[0] = 0;      /* X component */
    v2[1] = size;   /* Y component */
    v2[2] = c-b;    /* Z component */

    /* Obtain cross-product */
    crossp(v1,v2,norm);

    return;
} 



static void find_norm_lr(a,b,c,size,norm)
float a,b,c,size,*norm;
{
    /* Compute normal for lower right corner geometry */
    float v1[3],v2[3];

    /* Load the vectors */
    v1[0] = 0;      /* X component */
    v1[1] = -size;   /* Y component */
    v1[2] = a-b;    /* Z component */
    v2[0] = -size;   /* X component */
    v2[1] = 0;      /* Y component */
    v2[2] = c-b;    /* Z component */

    /* Obtain cross-product */
    crossp(v1,v2,norm);

    return;
}



static void find_norm_ll(a,b,c,size,norm)
float a,b,c,size,*norm;
{
    /* Compute normal for lower left corner geometry */
    float v1[3],v2[3];

    /* Load the vectors */
    v1[0] = size;   /* X component */
    v1[1] = 0;      /* Y component */
    v1[2] = a-b;    /* Z component */
    v2[0] = 0;      /* X component */
    v2[1] = -size;  /* Y component */
    v2[2] = c-b;    /* Z component */

    /* Obtain cross-product */
    crossp(v1,v2,norm);

    return;
}



#define MAG(x0,x1,x2) sqrt(x0*x0+x1*x1+x2*x2) 

static void crossp(v1,v2,norm)
float *v1,*v2,*norm;
{
    float m;
    /* Compute cross product in a LEFT HANDED COORDINATE SYSTEM */
    /* Form cross product norm = v1 X v2 */
    norm[0] =  -v1[1]*v2[2]+v1[2]*v2[1];
    norm[1] =  -v1[0]*v2[2]+v2[0]*v1[2];
    norm[2] =  -v1[0]*v2[1]+v2[0]*v1[1];
    m = sqrt((double)(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]));
    norm[0] /= m;
    norm[1] /= m;
    norm[2] /= m;

    return;
}
/* -library_code_end */
