Explore the World of Soft…e: Engineering & Science/Explore_the_World_of_Software_Engineering_and_Science_HRS_Software_1998.iso/programs/civil_en/dome350.exe/DOME.CPP

C/C++ Source or Header | 1995-04-08 | 47KB | 1,466 lines

//$Id: dome.cpp 3.50 1995/04/08 22:16:05 PONDERMATI Released PONDERMATI $

/*
DOME.CPP - A program & C++ class for the calculation of geodesic dome
            properties

     Copyright (C) 1995  Richard J. Bono

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

     Please direct inquiries, comments and modifications to:
     Richard J. Bono
     44 Augusta Rd.
     Brownsville, TX 78521

     email: [e-mail]
                -or-
            [e-mail]

Revision history:

$Log: dome.cpp $
'Revision 3.50  1995/04/08  22:16:05  PONDERMATI
'Added Class II symmetry triangle support
'Added spherical rotations of Class I polyhdera
'Added headers & spherical coordinates to PRN output
'Added routine to test bottom face tetra rotation
'Added rounding checks to rotate functions
'Class constructor passed a input parameter structure
'Removed interactive interface
'Streamlined verbose display
'Added separate help function
'Added Class I sphere formation
'Updated exit codes to non-zero where required
'Added non-Borland C defs
'Added automatic revision tracking
'Enabled Class II symmetry triangle support
'
'Revision 2.18  1995/02/12  04:54:49  PONDERMATIC
'Release Version of Dome.
'Added command line-support
'Added support for PRN & POV-Ray output files
'Rearranged many functions and parameter passing
'
'Revision 2.0  1995/01/01  22:54:59  PONDERMATI
'First public release. Supports only primary symmetry triangle calculations,
'class I (Alternate) breakdowns, & DXF & ASCII file outputs.
'

Acknowledgements & References:
The main reference used in the creation of this code was "Geodesic Math & How
to Use It" by Hugh Kenner, 1976, University of California Press.
ISBN 0-520-02924-0; Library of Congress Catalog Card Number: 74-27292. Many
thanks to Hugh for putting this data in an accessible format.

Also, many thanks to:
    J. F. Nystrom
    My wife Sonia
    My daughter Kathy
    Chris Fearnley
    Kirby Urner
        &
    R. Buckminster Fuller for changing the way I view Universe.

Compliation:
This program was complied and tested using Borland C++ 4.5 using the large
memory model on a Gateway 2000 4DX2-50V.

The release executables were complied using Borland C++ 3.1. The program will
now run on any IBM PC class machine with or withou a math coprocessor.

Support is included for non-Borland compilation. Unix patches provided thanks
to Chris Fearnley & John Kirk.
*/

#include<iostream.h>
#include<iomanip.h>
#include<fstream.h>
#include<stdlib.h>
#include<math.h>
#include<ctype.h>
#include<string.h>

#ifdef __BORLANDC__
#include<conio.h>
#endif

#ifdef __MSDOS__
#define HUGEPREFIX huge
#else
#define HUGEPREFIX
#endif

static char rcsid[]="$Id: dome.cpp 3.50 1995/04/08 22:16:05 PONDERMATI Released PONDERMATI $";

const double RAD_TO_DEG = 57.2957795130824;
const double DEG_TO_RAD = .0174532925199433;

#define BYTE unsigned char
#define WORD unsigned int

//Structure containing input parameters
struct parameters{
    long freq;
    long classt;
    long polyt;
    long filet;
    char filename[80];
    int verbose_flag;
    int sphere_flag;
} cmd_parm;

//Function prototypes
void logo_display(void);
void help_display(void);
void get_cmd(struct parameters *command, int param_count, char *param[]);
int main(int argc, char *argv[]);

//--------------geodesic points class
class Geodesic    {
private:
    //Various indices, and temporary variables
    long i,    j, k;
    long index;
    double sphi, stheta, ephi, etheta;
    //actual values used in calculations set dependent on class type
    long freq_calc, vertex_calc, face_calc, edges_calc;

public:
    long frequency, vertex, faces, edges;
    long sphere_flg, poly_face;
    long vertexII, facesII, edgesII;        //class II topology
    long classtype, polytype;

    //structure defines geodesic coordinate points
    struct coordinates{
        //vertex labels: 0,0 equals zenith
        long A;
        long B;
        //vertex coordinates: 0,0,0 equals zenith
        long xprime;
        long yprime;
        long zprime;

        //vertex spherical coordinates
        double phi;
        double theta;
        //Vertex cartesian coordinates
        double X;
        double Y;
        double Z;
    };

    //structure defines chord data
    struct chords{
        //start point
        long sA;
        long sB;
        double sphi;
        double stheta;

        //end point
        long eA;
        long eB;
        double ephi;
        double etheta;
    };

    //portability note: huge keyword was needed in order to allocate memory for
    //arrays from the far heap in a MSDOS far memory model.
    coordinates HUGEPREFIX *pntcrd;  //pointer to coordinate array
    chords HUGEPREFIX *edgepts;      //pointer to chord array

    Geodesic(struct parameters *command);    //constructor
    ~Geodesic();                               //destructor
    void topology(void);                       //topological abundance function
    void init_crd(void);                       //Set up coordinate values
    void spherical(void);                       //calculate spherical coordinates
    double chord_length(double, double, double, double); //chord length function
    void chord_factor(void);                   //calculates chord factors
    void save_dxf(char *);                     //save chord data in DXF format
    void save_ascii(char *);                   //save all coordinate data in ASCII
    void save_POV(char *);                    //save data in POV format
    void save_PRN(char *);                    //Save raw data in ASCII
    void display_data(void);                //Display data during program execution
    double clean_float(double);                //function cleans up triag zeros
    double rotate_phi(double, double, double, double);    //phi rotation function
    double rotate_theta(double, double, double, double); //theta rotation function
    void make_sphere(void);                    //generate data for a full sphere
    long tetra_angle(void);            //Get "A" coordinate to begin correction of bottom tetra face
};

//--------------geodesic constructor
Geodesic::Geodesic(struct parameters *command)
{

    //The constructor handles all calculation when an instance is called.
    //The input parameters are the frequency of subdivision, the class
    //type, the polygon type and a flag to determine whether to generate
    //spherical data. This is all that is needed tospecify a geodesic dome!
    //(Note: I have not included radius as size is special case and scaling
    //factors can be applied. The program assumes unit radius domes).

    //initialize class variables
    frequency = command->freq;
    classtype = command->classt;
    polytype = command->polyt;
    sphere_flg = command->sphere_flag;
    vertex = 0;
    faces = 0;
    edges = 0;

    //start calculations
    topology();            //calculate symmetry triangle topological abundance
    init_crd();            //initialize vertex coordinates and data structures
    spherical();         //calculate spherical coordinates
    chord_factor();      //determine chord factors
    if(sphere_flg)
        make_sphere();    //Generate full sphere coordinates if required

}
//--------------geodesic destructor
Geodesic::~Geodesic() {}
//-------------------calculate chord factors
void Geodesic::chord_factor(void)
{
    //Function cycles through each vertex and determines the
    //chord connections. This code is not very elegant but I think
    //it provides more readablity. Future incarnations may modify this to
    //reduce source code size.

    index = 1;
    for(i=1; i<vertex_calc; i++){
        if(pntcrd[i].A == 0 && pntcrd[i].B == 0){
            //point is the zenith vertex; add 10 & 11
            //line #1 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

            //line #2 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

        }
        else if(pntcrd[i].A == freq_calc && pntcrd[i].B == 0){
            //point is a left vertex; add 01
            //line #1 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;
        }
        else if(pntcrd[i].A == pntcrd[i].B){
            //point is a right edge; add 10 & 11
            //line #1 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

            //line #2 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;
        }
        else if(pntcrd[i].A == freq_calc){
            //point is a bottom vertex; add 01
            //line #1 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A  && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;
        }
        else if(pntcrd[i].B == 0){
            //point is a right edge; add 01, 10 & 11
            //line #1 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

            //line #2 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

            //line #3 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;
        }
        else {
            //point is an interior vertex; add 01, 10 & 11
            //line #1 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

            //line #2 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

            //line #3 start point definition
            edgepts[index].sA = pntcrd[i].A;
            edgepts[index].sB = pntcrd[i].B;
            edgepts[index].sphi = pntcrd[i].phi;
            edgepts[index].stheta = pntcrd[i].theta;

            for(j=i; j<=vertex_calc; j++){
                if(pntcrd[j].A == pntcrd[i].A + 1 && pntcrd[j].B == pntcrd[i].B + 1)
                    break;
            }
            //end point
            edgepts[index].eA = pntcrd[j].A;
            edgepts[index].eB = pntcrd[j].B;
            edgepts[index].ephi = pntcrd[j].phi;
            edgepts[index].etheta = pntcrd[j].theta;

            index++;

        }
    }

}
//--------------------clean floating point numbers
double Geodesic::clean_float(double fp_number)
{
        //function cleans floating point results close to but not quite zero
        if((fp_number < 0.0) & (fp_number > -0.0000000000001))
            fp_number = fabs(fp_number);

        return fp_number;
}
//-------------------save ASCII PRN format
void Geodesic::save_PRN(char *filename)
{
    double sX, sY, sZ, eX, eY, eZ;

    //This function saves the spherical & cartesian data to an
    //comma-delimited ASCII PRN file. This file can be imported into
    //most spreadsheets an the like.
    ofstream PRN(filename);

    //Set field widths
    PRN << setiosflags(ios::fixed) << setw(8) << setprecision(6);

    PRN << "Spherical Vertexia Coordinates" << '\n';
    //output PRN data...start with vertexia in spherical coordinates...
    for(j=1; j <= poly_face; j++){
        PRN << "Polyhedron Face #" << j << '\n';
        for(i=(vertex_calc * (j-1)) + 1; i <= (vertex_calc * j); i++)
            //Save data
            PRN << pntcrd[i].phi << ", " << pntcrd[i].theta << '\n';
    }
    PRN << '\n';

    PRN << "Cartesian Vertexia Coordinates" << '\n';
    //...then in XYZ coordinates
    for(j=1; j <= poly_face; j++){
        PRN << "Polyhedron Face #" << j << '\n';
        for(i=(vertex_calc * (j-1)) + 1; i <= (vertex_calc * j); i++){
            //convert spherical to cartesian
            sX = clean_float(cos(pntcrd[i].phi * DEG_TO_RAD) * sin(pntcrd[i].theta * DEG_TO_RAD));
            sY = clean_float(sin(pntcrd[i].phi * DEG_TO_RAD) * sin(pntcrd[i].theta * DEG_TO_RAD));
            sZ = clean_float(cos(pntcrd[i].theta * DEG_TO_RAD));
            //Save data
            PRN << sX << ", " << sY << ", " << sZ << '\n';
        }
    }

    PRN << '\n';

    PRN << "Spherical Chord Coordinates" << '\n';
    //Now do the chords...first in spherical coordinate pairs...
    for(j=1; j <= poly_face; j++){
        PRN << "Polyhedron Face #" << j << '\n';
        for(i=(edges_calc * (j-1)) + 1; i <= (edges_calc * j); i++)
            //save data
            PRN << edgepts[i].sphi << ", " << edgepts[i].stheta << ", "
                << edgepts[i].ephi << ", " << edgepts[i].etheta << '\n';
    }

    PRN << '\n';

    PRN << "Cartesian Chord Coordinates" << '\n';
    //...then in XYZ...
    for(j=1; j <= poly_face; j++){
        PRN << "Polyhedron Face #" << j << '\n';
        for(i=(edges_calc * (j-1)) + 1; i <= (edges_calc * j); i++){
            //convert spherical to cartesian
            sX = clean_float(cos(edgepts[i].sphi * DEG_TO_RAD) * sin(edgepts[i].stheta * DEG_TO_RAD));
            sY = clean_float(sin(edgepts[i].sphi * DEG_TO_RAD) * sin(edgepts[i].stheta * DEG_TO_RAD));
            sZ = clean_float(cos(edgepts[i].stheta * DEG_TO_RAD));
            eX = clean_float(cos(edgepts[i].ephi * DEG_TO_RAD) * sin(edgepts[i].etheta * DEG_TO_RAD));
            eY = clean_float(sin(edgepts[i].ephi * DEG_TO_RAD) * sin(edgepts[i].etheta * DEG_TO_RAD));
            eZ = clean_float(cos(edgepts[i].etheta * DEG_TO_RAD));
            //save data
            PRN << sX << ", " << sY << ", " << sZ << ", "
                 << eX << ", " << eY << ", " << eZ << '\n';
        }
    }

    PRN.close();

}
//-------------------save POV-Ray file
void Geodesic::save_POV(char *filename)
{
    double sX, sY, sZ, eX, eY, eZ;

    //save chord data for symmetry triangle to POV file
    //filename should include .POV extension on MS-DOS systems

    ofstream POV(filename);
    //Set field widths
    POV << setiosflags(ios::fixed) << setw(8) << setprecision(6);

    //set up file header
    POV << "//POV-Ray script" << '\n';
    POV << "#include" << '"' << "colors.inc" << '"' << '\n' << '\n';
    POV << "#declare Cam_factor = 6" << '\n';
    POV << "#declare Camera_X = 1 * Cam_factor" << '\n';
    POV << "#declare Camera_Y = 0.5 * Cam_factor" << '\n';
    POV << "#declare Camera_Z = -0.3 * Cam_factor" << '\n';

    POV << "camera { location  <Camera_X, Camera_Y, Camera_Z>" << '\n';
    POV << "        up        <0, 1.0,  0>    right     <-1.33, 0,  0>" << '\n';
    POV << "        direction <0, 0,  3>      look_at   <0, 0, 0> }" << '\n' << '\n';

    POV << "light_source { <Camera_X - 2, Camera_Y + 5 , Camera_Z + 5> color White }" << '\n';
    POV << "light_source { <Camera_X - 2, Camera_Y + 5 , Camera_Z - 3> color White }" << '\n';

    POV << '\n' << "// Background:" << '\n';
    POV << "background {color Gray70}" << '\n' << '\n';

    POV << "// POV script for DOME: " << '\n' << '\n';

    POV << "union {" << '\n';

    //output POV data...start with spheres at vertexia points...
    for(i=1; i<=vertex_calc * poly_face; i++){
        //convert spherical to cartesian
        sX = clean_float(cos(pntcrd[i].phi * DEG_TO_RAD) * sin(pntcrd[i].theta * DEG_TO_RAD));
        sY = clean_float(sin(pntcrd[i].phi * DEG_TO_RAD) * sin(pntcrd[i].theta * DEG_TO_RAD));
        sZ = clean_float(cos(pntcrd[i].theta * DEG_TO_RAD));
        //Save data
        POV << "sphere{<" << sX << ", " << sY << ", " << sZ;
        POV << ">,0.04 pigment {Blue} no_shadow}" << '\n';
    }
    POV << '\n';

    //Now do the chords
    for(i=1; i<=edges_calc * poly_face; i++){
        //convert spherical to cartesian
        sX = clean_float(cos(edgepts[i].sphi * DEG_TO_RAD) * sin(edgepts[i].stheta * DEG_TO_RAD));
        sY = clean_float(sin(edgepts[i].sphi * DEG_TO_RAD) * sin(edgepts[i].stheta * DEG_TO_RAD));
        sZ = clean_float(cos(edgepts[i].stheta * DEG_TO_RAD));
        eX = clean_float(cos(edgepts[i].ephi * DEG_TO_RAD) * sin(edgepts[i].etheta * DEG_TO_RAD));
        eY = clean_float(sin(edgepts[i].ephi * DEG_TO_RAD) * sin(edgepts[i].etheta * DEG_TO_RAD));
        eZ = clean_float(cos(edgepts[i].etheta * DEG_TO_RAD));
        //save data
        POV << "cylinder{<" << sX << ", " << sY << ", " << sZ << ">, ";
        POV << "<" << eX << ", " << eY << ", " << eZ << ">,0.03 pigment {Red} no_shadow}" << '\n';
    }
    POV << "rotate <0, 0, 0> }" << '\n';

    POV.close();

}
//-------------------save DXF file
void Geodesic::save_dxf(char *filename)
{
    double sX, sY, sZ, eX, eY, eZ;

    //save chord data for symmetry triangle to DXF file
    //filename should include .DXF extension on MS-DOS systems
    //For full spheres each face is saved on a different level...

    ofstream DXF(filename);
    //output DXF data
    DXF << "0" << '\n';
    DXF << "SECTION" << '\n';
    DXF << "2" << '\n';
    DXF << "ENTITIES" << '\n';
    //Set field widths
    DXF << setiosflags(ios::fixed) << setw(8) << setprecision(6);

    for(j=1; j<=poly_face; j++){
        for(i=(edges_calc * (j-1)) + 1; i<=(edges_calc * j); i++){
            //convert spherical to cartesian
            sX = clean_float(cos(edgepts[i].sphi * DEG_TO_RAD) * sin(edgepts[i].stheta * DEG_TO_RAD));
            sY = clean_float(sin(edgepts[i].sphi * DEG_TO_RAD) * sin(edgepts[i].stheta * DEG_TO_RAD));
            sZ = clean_float(cos(edgepts[i].stheta * DEG_TO_RAD));
            eX = clean_float(cos(edgepts[i].ephi * DEG_TO_RAD) * sin(edgepts[i].etheta * DEG_TO_RAD));
            eY = clean_float(sin(edgepts[i].ephi * DEG_TO_RAD) * sin(edgepts[i].etheta * DEG_TO_RAD));
            eZ = clean_float(cos(edgepts[i].etheta * DEG_TO_RAD));
            //save data
            DXF << "0" << '\n';
            DXF << "LINE" << '\n';
            DXF << "8" << '\n';
            DXF << j << '\n';
            DXF << "10" << '\n';
            DXF << sX << '\n';
            DXF << "20" << '\n';
            DXF << sY << '\n';
            DXF << "30" << '\n';
            DXF << sZ << '\n';
            DXF << "11" << '\n';
            DXF << eX << '\n';
            DXF << "21" << '\n';
            DXF << eY << '\n';
            DXF << "31" << '\n';
            DXF << eZ << '\n';
        }
    }
    DXF << "0" << '\n';
    DXF << "ENDSEC" << '\n';
    DXF << "0" << '\n';
    DXF << "EOF" << '\n';

    DXF.close();

}
//-------------------save data to ASCII file
void Geodesic::save_ascii(char *filename)
{
    //function to save data to ASCII data file
    //WARNING: CAN produce very large listings!
    //save chord data for symmetry triangle to DAT file
    //filename should include .DAT extension on MS-DOS systems
    //Only data for the symmetry triangle is saved

    double sphi, stheta, ephi, etheta;

    ofstream ASCII(filename);
    //output ASCII report data
    ASCII << "Geodesic Dome Data" << '\n';
    ASCII << "----------------Dome Parameters-------------------------" << '\n';
    ASCII << "Polyhedron Type: ";
    if(polytype == 1)
        ASCII << "Icosahedron" << '\n';
    else if(polytype == 2)
        ASCII << "Octahedron" << '\n';
    else
        ASCII << "Tetrahedron" << '\n';
    ASCII << "Class ";
    if(classtype == 1)
        ASCII << "I" << '\n';
    else
        ASCII << "II" << '\n';
    ASCII << "Frequency: " << frequency << '\n';
    ASCII << "----------------Symmetry Triangle Data------------------" << '\n';
    ASCII << "Vertexia: " << vertex << '\n';
    ASCII << "Edges: " << edges << '\n';
    ASCII << "Faces: " << faces << '\n';
    ASCII << "-----------Symmetry Triangle Spherical Coordinates------" << '\n';
    ASCII << "A" << '\t' << "B" << '\t' << "phi" << '\t' << "theta" << '\n';

    for(i=1; i<=vertex_calc; i++){
        ASCII << setiosflags(ios::right) << setw(3);
        ASCII << pntcrd[i].A << '\t';
        ASCII << pntcrd[i].B << '\t';
        ASCII << setw(11) << setprecision(8) << setiosflags(ios::fixed);
        ASCII << pntcrd[i].phi << '\t';
        ASCII << pntcrd[i].theta << '\n';
    }
    ASCII << '\n';
    ASCII << "----------Symmetry Triangle Chord Data------------------" << '\n';
    ASCII << "AB start/AB end" << '\t' << "Chord Length" << '\n';
    for(i=1; i<=edges_calc; i++){
        ASCII << setiosflags(ios::right) << setw(3);
        ASCII << edgepts[i].sA << "," << edgepts[i].sB << "/";
        ASCII << edgepts[i].eA << "," << edgepts[i].eB << '\t';
        ASCII << setw(12) << setprecision(8) << setiosflags(ios::fixed);
        sphi = edgepts[i].sphi;
        stheta = edgepts[i].stheta;
        ephi = edgepts[i].ephi;
        etheta = edgepts[i].etheta;
        ASCII << chord_length(sphi, stheta, ephi, etheta) << '\n';
    }
    ASCII << "End Dome Data" << '\n';

    ASCII.close();
}
//-------------------calculate chord lengths
double Geodesic::chord_length(double sphi, double stheta, double ephi, double etheta)
{
    //Function returns distance between two points given their
    //spherical coordinates. Radius is fixed at one.

    double result;

    result = pow((2-2 * (cos(stheta * DEG_TO_RAD) * cos(etheta * DEG_TO_RAD) +
        cos((sphi-ephi) *DEG_TO_RAD) * sin(stheta * DEG_TO_RAD) *
            sin(etheta * DEG_TO_RAD))),0.5);

    return result;

}
//------------------phi rotation function
double Geodesic::rotate_phi(double phi, double phi1, double theta1, double theta2)
{
        //phi half of the rotation formula given in Kenner's "Geodesic Math &
        //How to Use it" Kenner credits this formula to Professor George Owen

        double result;

        result = (sin(theta1 * DEG_TO_RAD)* sin((phi-phi1) * DEG_TO_RAD)) /
                        sin(theta2 * DEG_TO_RAD);

        //Apply correction for round-off errors
        if(result > 1.0)
            result = 1.0;
        else if(result < -1.0)
            result = -1.0;

        return asin(result) * RAD_TO_DEG;

}
//------------------theta rotation function
double Geodesic::rotate_theta(double phi, double phi1, double theta, double theta1)
{
        //theta half of the rotation formula given in Kenner's "Geodesic Math &
        //How to Use it" Kenner credits this formula to Professor George Owen
        double result;

        result = cos(theta1 * DEG_TO_RAD) * cos(theta * DEG_TO_RAD) +
                        sin(theta1 * DEG_TO_RAD) * sin(theta * DEG_TO_RAD) *
                            cos((phi-phi1) * DEG_TO_RAD);

        //Apply correction for round-off errors
        if(result > 1.0)
            result = 1.0;
        else if(result < -1.0)
            result = -1.0;

        return acos(result) * RAD_TO_DEG;

}
//--------------calculate topological abundance of symmetry triangle
void Geodesic::topology(void)
{

    //This class function calculates the topological abundance
    //of the symmetry triangle.

    //calculate the number of vertexia & faces in symmetry triangle
    //This is based on the class type
    //First the polyhedron face...

    for(i=1; i<=frequency; i++){
        if(i == 1)
            vertex = 3;
        else
            vertex += (i + 1);
    }
    faces = frequency * frequency;
    edges = (vertex + faces) - 1;

    //if class II then calculate the class II symmetry triangle topology.
    //The class II triangle is NOT the polyhedron face!
    if(classtype == 2){
        for(i=1; i<=frequency/2; i++){
            if(i == 1)
                vertexII = 3;
            else
                vertexII += (i + 1);
        }
        facesII = (frequency/2) * (frequency/2);
        edgesII = (vertexII + facesII) - 1;
    }
}
//-------------initialize coordinate array by setting vertex points
void Geodesic::init_crd(void)
{
    //this class function first allocates memory for the coordinate
    //and chord arrays. The vertex labels and vertex coordinates are then
    //calculated. See chapter 12 of Kenner's "Geodesic Math & How to Use
    //it" for more specifics on the labeling & coordinates

    //Set the topology dependent on class type
    if(classtype == 1){
        vertex_calc = vertex;
        edges_calc = edges;
        freq_calc = frequency;
    }
    else if(classtype == 2){
        vertex_calc =vertexII;
        edges_calc = edgesII;
        freq_calc = frequency / 2;
    }


    //Determine number of faces required if a full sphere is to be generated
    if(sphere_flg){
        if(polytype == 1)
            poly_face = 20;        //Icosa = 20 faces
        else if(polytype == 2)
            poly_face = 8;        //Octa = 8 faces
        else
            poly_face = 4;        //Tetra = 4 faces
    }
    else
        poly_face = 1;            //Symmetry triangle only

    //Memory is allocated for the full polyhedron face. Create array object for vertexia
    //In the case of a sphere, memory is allocated for everything
    pntcrd = new HUGEPREFIX coordinates[(vertex + 1) * poly_face];
    if (pntcrd == NULL){
        logo_display();
        cout << "Insufficient memory for coordinate array --- Execution Terminating" << '\n';
        exit(1);
    }

    //create array object for edges
    edgepts = new HUGEPREFIX chords[(edges + 1) * poly_face];
    if (edgepts == NULL){
        logo_display();
        cout << "Insufficient memory for chord array --- Execution Terminating" << '\n';
        exit(1);
    }

    //Initialize vertex labels and coordinates
    index = 1;
    for(i=0; i<=frequency; i++){
        for(j=0; j<=i; j++){
            pntcrd[index].A = i;
            pntcrd[index].B = j;
            pntcrd[index].xprime = j;
            pntcrd[index].yprime = i - j;
            pntcrd[index].zprime = frequency - i;
            index++;
        }
    }
}
//-------------Determine "A" Coordinate of bottom tetra face-----------
long Geodesic::tetra_angle(void)
{
        //This function returns the "A" coordinate of the bottom tetrahedron face
        //where correction of the phi angle must begin. It is the result of quirks
        //in professor owens rotation formulae. I don't find this code very elegant
        //but it does provide the proper results.

        double sX, sphi, stheta;
        double eX, ephi, etheta;

        //Save the 0,0 coordinate as the start point and convert to an X value
        stheta = rotate_theta(240.0, pntcrd[1].phi, acos(-1.0/3.0) * RAD_TO_DEG,
                        pntcrd[1].theta);
        sphi = rotate_phi(240.0, pntcrd[1].phi, pntcrd[1].theta, stheta);
        sX = clean_float(cos(sphi * DEG_TO_RAD) * sin(stheta * DEG_TO_RAD));

        //The goal of this routine is to find the "A" coordinate in which the
        //difference between the current X and the last X is positive. We start
        //at 0,0 then proceed down the tetra face edge until this is found
        //(i.e. 0,0 to 1,0; 1,0 to 2,0; 2,0 to 3,0;...f-1,0 to f,0)

        index = 1;
        for(i=2; i<=freq_calc; i++){
            index += (i-1);    //index given the array location of the tetra edge point

            //Save the next coordinate as the end point and convert to an X value
            etheta = rotate_theta(240.0, pntcrd[index].phi, acos(-1.0/3.0) * RAD_TO_DEG,
                            pntcrd[index].theta);
            ephi = rotate_phi(240.0, pntcrd[index].phi, pntcrd[index].theta, etheta);

            eX = clean_float(cos(ephi * DEG_TO_RAD) * sin(etheta * DEG_TO_RAD));

            if((eX - sX) > 0.0)
                return pntcrd[index].A;
            else
                sX = eX;
        }

        //Handle special cases of frequency < 5
        if(freq_calc == 4)
            return 3;
        if(freq_calc == 1)
            return 1;
        else
            return 2;

}
//-------------rotate symmetry triangle into sphere--------------------
void Geodesic::make_sphere(void)
{
    //Rotates chords & points to provide for full sphere
    //Set the topology dependent on class type
    if(classtype == 1){
        vertex_calc = vertex;
        edges_calc = edges;
    }
    else if(classtype == 2){
        vertex_calc =vertexII;
        edges_calc = edgesII;
    }

    if(polytype == 1 && classtype == 1){
        //Do downward pointing triangle
        index = edges_calc + 1;
        for(j=1; j<=edges_calc; j++){
            edgepts[index].stheta = rotate_theta(36.0, edgepts[j].sphi, (180.0 - (atan(2) * RAD_TO_DEG)),
                        edgepts[j].stheta);
            edgepts[index].sphi = (rotate_phi(36.0, edgepts[j].sphi, edgepts[j].stheta,
                        edgepts[index].stheta) + 36.0);
            edgepts[index].etheta = rotate_theta(36.0, edgepts[j].ephi, (180.0 - (atan(2) * RAD_TO_DEG)),
                        edgepts[j].etheta);
            edgepts[index].ephi = (rotate_phi(36.0, edgepts[j].ephi, edgepts[j].etheta,
                        edgepts[index].etheta) + 36.0);
            index++;
        }
        //Rotate triangles #1 & #2 4 times
        for(i=1; i<=4; i++){
            for(j=1; j<=edges_calc * 2; j++){
                  edgepts[index].sphi = edgepts[j].sphi + 72.0 * i;
                  edgepts[index].stheta = edgepts[j].stheta;
                  edgepts[index].ephi = edgepts[j].ephi + 72.0 * i;
                  edgepts[index].etheta = edgepts[j].etheta;
                  index++;
            }
        }
        //Rotate 10 triangles by 180 degrees of theta
        for(j=1; j<=edges_calc * 10; j++){
              edgepts[index].sphi = edgepts[j].sphi + 36.0;
              edgepts[index].stheta = 180.0 - edgepts[j].stheta;
              edgepts[index].ephi = edgepts[j].ephi + 36.0;
              edgepts[index].etheta = 180.0 - edgepts[j].etheta;
              index++;
        }

        //now repeat for the vertexia
        //Do downward pointing triangle
        index = vertex_calc + 1;
        for(j=1; j<=vertex_calc; j++){
            pntcrd[index].theta = rotate_theta(36.0, pntcrd[j].phi, (180.0 - (atan(2) * RAD_TO_DEG)),
                        pntcrd[j].theta);
            pntcrd[index].phi = (rotate_phi(36.0, pntcrd[j].phi, pntcrd[j].theta,
                        pntcrd[index].theta) + 36.0);
            index++;
        }
        //Rotate triangles #1 & #2 4 times
        for(i=1; i<=4; i++){
            for(j=1; j<=vertex_calc * 2; j++){
                  pntcrd[index].phi = pntcrd[j].phi + 72.0 * i;
                  pntcrd[index].theta = pntcrd[j].theta;
                  index++;
            }
        }
        //Rotate 10 triangles by 180 degrees of theta
        for(j=1; j<=vertex_calc * 10; j++){
              pntcrd[index].phi = pntcrd[j].phi + 36.0;
              pntcrd[index].theta = 180.0 - pntcrd[j].theta;
              index++;
        }

    }
    else if(polytype == 2 && classtype == 1){
        //Rotate Octahedron; Class II
        //rotate chords to form hemisphere
        index = edges_calc + 1;
        for(i=1; i<=3; i++){
            for(j=1; j<=edges_calc; j++){
                  edgepts[index].sphi = edgepts[j].sphi + 90.0 * i;
                  edgepts[index].stheta = edgepts[j].stheta;
                  edgepts[index].ephi = edgepts[j].ephi + 90.0 * i;
                  edgepts[index].etheta = edgepts[j].etheta;
                  index++;
            }
        }
        //do next hemisphere
        for(i=0; i<=3; i++){
            for(j=1; j<=(edges_calc); j++){
                  edgepts[index].sphi = edgepts[j].sphi + 90.0 * i;
                  edgepts[index].stheta = 180.0 - edgepts[j].stheta;
                  edgepts[index].ephi = edgepts[j].ephi + 90.0 * i;
                  edgepts[index].etheta = 180.0 - edgepts[j].etheta;
                  index++;
            }
        }
        //rotate vertexia to form hemisphere
        index = vertex_calc + 1;
        for(i=1; i<=3; i++){
            for(j=1; j<=vertex_calc; j++){
                  pntcrd[index].phi = pntcrd[j].phi + 90.0 * i;
                  pntcrd[index].theta = pntcrd[j].theta;
                  index++;
            }
        }
        //do next hemisphere
        for(i=0; i<=3; i++){
            for(j=1; j<=(vertex_calc); j++){
                  pntcrd[index].phi = pntcrd[j].phi + 90.0 * i;
                  pntcrd[index].theta = 180.0 - pntcrd[j].theta;
                  index++;
            }
        }

    }
    else if(polytype == 3 && classtype == 1){
        //Rotate tetrahedron; Class I
        //find A coordinate where angles turn
        long A_change = tetra_angle();

        //rotate chords to form first three faces
        index = edges_calc + 1;
        for(i=1; i<=2; i++){
            for(j=1; j<=edges_calc; j++){
                  edgepts[index].sphi = edgepts[j].sphi + 120.0 * i;
                  edgepts[index].stheta = edgepts[j].stheta;
                  edgepts[index].ephi = edgepts[j].ephi + 120.0 * i;
                  edgepts[index].etheta = edgepts[j].etheta;
                  index++;
            }
        }
        //use rotation formula for bottom face chords
        for(j=1; j<=edges_calc; j++){
            edgepts[index].stheta = rotate_theta(240.0, edgepts[j].sphi, acos(-1.0/3.0) * RAD_TO_DEG,
                        edgepts[j].stheta);
            edgepts[index].sphi = rotate_phi(240.0, edgepts[j].sphi, edgepts[j].stheta,
                        edgepts[index].stheta);
            edgepts[index].etheta = rotate_theta(240.0, edgepts[j].ephi, acos(-1.0/3.0) * RAD_TO_DEG,
                        edgepts[j].etheta);
            edgepts[index].ephi = rotate_phi(240.0, edgepts[j].ephi, edgepts[j].etheta,
                        edgepts[index].etheta);

            if(edgepts[j].sA >= A_change)
                edgepts[index].sphi = 180.0 - edgepts[index].sphi;
            if(edgepts[index].sphi < 0.0)
                edgepts[index].sphi += 360.0;
            if(edgepts[j].eA >= A_change)
                edgepts[index].ephi = 180.0 - edgepts[index].ephi;
            if(edgepts[index].ephi < 0.0)
                edgepts[index].ephi += 360.0;

            index++;
        }

        //now do the vertexia...
        index = vertex_calc + 1;
        for(i=1; i<=2; i++){
            for(j=1; j<=vertex_calc; j++){
                  pntcrd[index].phi = pntcrd[j].phi + 120.0 * i;
                  pntcrd[index].theta = pntcrd[j].theta;
                  index++;
            }
        }
        //use rotation formula for bottom face vertexia
        for(j=1; j<=vertex_calc; j++){
            pntcrd[index].theta = rotate_theta(240.0, pntcrd[j].phi, acos(-1.0/3.0) * RAD_TO_DEG,
                        pntcrd[j].theta);
            pntcrd[index].phi = rotate_phi(240.0, pntcrd[j].phi, pntcrd[j].theta,
                        pntcrd[index].theta);
            if(pntcrd[j].A >= A_change)
                pntcrd[index].phi = 180.0 - pntcrd[index].phi;
            if(pntcrd[index].phi < 0.0)
                pntcrd[index].phi += 360.0;

            index++;
        }
    }
    else if(polytype == 1 && classtype == 2){
        //Rotate Icosahedra; Class II
        logo_display();
        cout << "Class II Sphere Not Supported --- Execution Terminating" << '\n';
        exit(1);
    }
    else if(polytype == 2 && classtype == 2){
        //Rotate ocatahedron; Class II
        logo_display();
        cout << "Class II Sphere Not Supported --- Execution Terminating" << '\n';
        exit(1);
    }
    else if(polytype == 3 && classtype == 2){
        //Rotate tetrahedron; Class II
        logo_display();
        cout << "Class II Sphere Not Supported --- Execution Terminating" << '\n';
        exit(1);
    }
}
//-------------calculate spherical coordinates for given polygon type & class
void Geodesic::spherical(void)
{
    //calculate spherical coordinates of geodesic vertexia
    //given coordinates values
    //care must be taken to avoid generating a trig error

    //Set the topology dependent on class type
    if(classtype == 1)
        vertex_calc = vertex;
    else if(classtype == 2)
        vertex_calc =vertexII;

    for(i=1; i<=vertex_calc; i++){
        //set-up x,y,z coordinates for spherical calulations
        //See page 75 of Kenner's "Geodesic Math & How to use it"
        //for formulae used below
        if(polytype == 2){
                //Octahedron
                pntcrd[i].X = pntcrd[i].xprime;
                pntcrd[i].Y = pntcrd[i].yprime;
                pntcrd[i].Z = pntcrd[i].zprime;
        }else if(polytype == 1){
                //Icosahedron
                pntcrd[i].X = pntcrd[i].xprime * sin(72.0 * DEG_TO_RAD);
                pntcrd[i].Y = pntcrd[i].yprime + (pntcrd[i].xprime * cos(72.0 * DEG_TO_RAD));
                pntcrd[i].Z = frequency/2.0 + (pntcrd[i].zprime / (2.0 * cos(36.0 * DEG_TO_RAD)));
        }else if(polytype == 3){
                //Tetrahedron
                pntcrd[i].X = pntcrd[i].xprime * pow(3.0, 0.5);
                pntcrd[i].Y = 2 * pntcrd[i].yprime - pntcrd[i].xprime;
                pntcrd[i].Z = (3.0 * pntcrd[i].zprime - pntcrd[i].xprime - pntcrd[i].yprime) / pow(2.0, 0.5);
        }


        //Calculate phi while avoiding trig errors
        if(pntcrd[i].Y == 0 && pntcrd[i].X == 0)
            pntcrd[i].phi = 0.0;
        else if(pntcrd[i].Y == 0)
            pntcrd[i].phi = 90.0;
        else
            pntcrd[i].phi = atan(pntcrd[i].X / pntcrd[i].Y) * RAD_TO_DEG;

        //adjust value to correct quadrant
        if(pntcrd[i].phi < 0.0)
            pntcrd[i].phi += 180.0;

        //now calculate theta...this is class dependent...
        if(pntcrd[i].Z == 0)
            pntcrd[i].theta = 90.0;
        else if(classtype == 1)
            //All class I types use this equation for theta
            pntcrd[i].theta = atan((pow(pow(pntcrd[i].X, 2.0) + pow(pntcrd[i].Y, 2.0), 0.5)) / pntcrd[i].Z) * RAD_TO_DEG;
        else if(polytype == 2 && classtype == 2)
            //theta for class II octahedra
            pntcrd[i].theta = atan((pow((2.0 * (pow(pntcrd[i].X, 2.0) + pow(pntcrd[i].Y, 2.0))), 0.5)) / pntcrd[i].Z) * RAD_TO_DEG;
        else if(polytype == 1 && classtype == 2)
            //theta for class II Icosahedron
            pntcrd[i].theta = atan((pow(pow(pntcrd[i].X, 2.0) + pow(pntcrd[i].Y, 2.0), 0.5)) / (cos(36.0 * DEG_TO_RAD) * pntcrd[i].Z)) * RAD_TO_DEG;
        else if(polytype == 3 && classtype ==2)
            //theta for class II tetrahedron
            //The formula given in Geodesic math (eq 12.15) appears to be inccorrect. The formula
            //below provides results indicated on tables
            pntcrd[i].theta = atan((2 * pow((pow(pntcrd[i].X, 2.0) + pow(pntcrd[i].Y, 2.0)), 0.5) / pntcrd[i].Z)) * RAD_TO_DEG;

        //make sure the right quadrant is used
        if(pntcrd[i].theta < 0.0)
            pntcrd[i].theta += 180.0;
    }
}
//--------------------display data function
void Geodesic::display_data(void)
{
        //Set the topology dependent on class type
        if(classtype == 1){
            vertex_calc = vertex;
            edges_calc = edges;
        }
        else if(classtype == 2){
            vertex_calc =vertexII;
            edges_calc = edgesII;
        }
        cout << "Geodesic Dome Data" << '\n';
        cout << "----------------Dome Parameters-------------------------" << '\n';
        cout << "Polyhedron Type: ";
        if(polytype == 1)
            cout << "Icosahedron" << '\n';
        else if(polytype == 2)
            cout << "Octahedron" << '\n';
        else
            cout << "Tetrahedron" << '\n';
        cout << "Class ";
        if(classtype == 1)
            cout << "I" << '\n';
        else
            cout << "II" << '\n';
        cout << "Frequency: " << frequency << '\n';
        cout << "----------------Symmetry Triangle Data------------------" << '\n';
        cout << "Vertexia: " << vertex << '\n';
        cout << "Edges: " << edges << '\n';
        cout << "Faces: " << faces << '\n' << '\n';

}
//--------------------End of class Geodesic

void help_display(void)
{
    //Display usage and help then exit
    logo_display();
    cout << "Usage: dome -fnnn -cn -px [-s] [-v] [-h] [filename.xxx]" << '\n';
    cout << "Where: -fnnn is geodesic frequency (default nnn=2)" << '\n';
    cout << "       -cn is class type (n=1 or 2; default n=1)" << '\n';
    cout << "       -px sets the polyhedron type" << '\n';
    cout << "        where x is: i for icosahedron (default)" << '\n';
    cout << "                    o for octahedron" << '\n';
    cout << "                    t for tetrahedron" << '\n';
    cout << "       -s generate full sphere data (default: symmetry triangle)" << '\n';
    cout << "       -v verbose data display at run-time" << '\n';
    cout << "       -h displays this help screen" << '\n';
    cout << "        filename.xxx is a standard DOS filename" << '\n';
    cout << "        where xxx is: DXF, DAT, POV or PRN" << '\n';
    exit(2);
}
void logo_display(void)
{
    char rev[]="$Revision: 3.50 $";
    int i, j;
    char level[6];

    //Get the revision level. I have been using GNU RCS to do revision control
    //on dome. This revision level will be used to automatically display the
    //correct rev level.
    j = 11;
    i = 0;
    while(rev[j] != ' '){
        level[i] = rev[j];
        j++;
        i++;
    }
    level[i] = '\0';

    #ifdef __BORLANDC__
    clrscr();
    #endif

    cout << "Dome " << level << ", Copyright (C) 1995, Richard J. Bono" << '\n';
    cout << "Dome comes with ABSOLUTELY NO WARRANTY. This is free software," << '\n';
    cout << "and you are welcome to redistribute it under certain conditions." << '\n';
    cout << "See GNU General Public License for more details." << '\n' << '\n';
}

void get_cmd(struct parameters *command, int param_count, char *param[])
{
    //Get and parse command line
    char cmd_parm[6];

    //Set defaults
    command->freq = 2;            //frequency = 2
    command->classt = 1;        //Class I
    command->polyt = 1;            //Icosahedron
    command->verbose_flag = 0;  //disable verbose
    command->filet = 5;            //No file output
    command->sphere_flag = 0;    //Symmetry triangle only

    int t, j, k;

    for(t=1; t<param_count; ++t){
        if(param[t][0] == '-'){
            switch (tolower(param[t][1])){
                case ('p'):
                    //Set polyhedron type
                    if(tolower(param[t][2]) == 'i')
                        //Set to Icosahedron
                        command->polyt = 1;
                    else if(tolower(param[t][2]) == 'o')
                        //Set to octahedron
                        command->polyt = 2;
                    else if(tolower(param[t][2]) == 't')
                        //Set to tetrahedron
                        command->polyt = 3;
                    else{
                        logo_display();
                        cout << "Invalid Polyhedron Type --- Execution Terminating" << '\n';
                        exit(1);
                    }
                    break;
                case ('h'):
                    //Display help and exit. This overrides other parameters
                    help_display();
                    break;
                case ('v'):
                    //Enable Parameter Display during execution
                    command->verbose_flag = 1;
                    break;
                case ('s'):
                    //generate sphere
                    command->sphere_flag = 1;
                    break;
                case ('f'):
                    //get frequency
                    j=2;
                    k=0;
                    while(param[t][j]){
                        cmd_parm[k] = param[t][j];
                        j++;
                        k++;
                    }
                    cmd_parm[k] = '\0';
                    command->freq = atol(cmd_parm);
                    if (command->freq <= 0){
                        logo_display();
                        cout << "Invalid Frequency --- Execution Terminating" << '\n';
                        exit(1);
                    }
                    break;
                case ('c'):
                    //get class type
                    j=2;
                    k=0;
                    while(param[t][j]){
                        cmd_parm[k] = param[t][j];
                        j++;
                        k++;
                    }
                    cmd_parm[k] = '\0';
                    command->classt = atol(cmd_parm);
                    if (command->classt < 1 || command->classt > 2){
                        logo_display();
                        cout << "Invalid Class Type --- Execution Terminating" << '\n';
                        exit(1);
                    }
                    break;
                default:
                    logo_display();
                    cout << "Invalid command-line --- Execution Terminating" << '\n';
                    exit(1);
                    break;
            }
        }
        else{
            //Check to see if parameter is a valid filename
            j = 0;
            while(param[t][j]){
                if(param[t][j] != '.')
                    j++;
                else{
                    //Check for valid extension
                    j++;
                    k = 0;
                    while(param[t][j]){
                        cmd_parm[k] = tolower(param[t][j]);
                        j++;
                        k++;
                    }
                    cmd_parm[k] = '\0';
                    if(!strcmp(cmd_parm,"dxf"))
                        command->filet = 1;
                    else if(!strcmp(cmd_parm, "dat"))
                        command->filet = 2;
                    else if(!strcmp(cmd_parm, "pov"))
                        command->filet = 3;
                    else if(!strcmp(cmd_parm, "prn"))
                        command->filet = 4;
                    else{
                        logo_display();
                        cout << "Invalid command-line --- Execution terminating" << '\n';
                        exit(1);
                    }
                    j = 0;
                    k = 0;
                    while(param[t][j]){
                        command->filename[k] = tolower(param[t][j]);
                        j++;
                        k++;
                    }
                    command->filename[k] = '\0';
                }
            }
        }
    }
    if(fmod(command->freq, 2) != 0 && command->classt == 2){
        logo_display();
        cout << "Class II requires Even Frequency --- Execution Terminating" << '\n';
        exit(1);
    }
}

int main(int argc, char *argv[])
{
    //This main routines shows what can be done with this class. It implements
    //a simple program which displays dome parameters and optionally saves the
    //data in either DXF or ASCII file formats.

    //Some things that can be tried involve creating several instances of the
    //geodesic object. Geodesic shells can be created in this way given enough
    //memory.

    //Command line parameters are changing this routine into something more
    //permanent. Remember the class can be incorporated into your own programs.
    //This shell is just one possible implementation.

    //Get dome parameters
    if(argc == 1)
        //display usage and exit
        help_display();
    else
        //Get command line
        get_cmd(&cmd_parm, argc, argv);

    //class instance
    Geodesic geosys(&cmd_parm);

    logo_display();
    if(cmd_parm.verbose_flag)
        geosys.display_data();

    if(cmd_parm.filet == 1)
        //output file in DXF format
        geosys.save_dxf(cmd_parm.filename);
    else if(cmd_parm.filet == 2)
        //output all data in ASCII report format
        geosys.save_ascii(cmd_parm.filename);
    else if(cmd_parm.filet ==3)
        //output data in POV-Ray script format
        geosys.save_POV(cmd_parm.filename);
    else if(cmd_parm.filet ==4)
        //output data in PRN format
        geosys.save_PRN(cmd_parm.filename);

    cout << '\n' << "Execution Complete!" << '\n';
    return(0);
}