//dawn.c  and szg-1.1 11jun07 by gkf 
//rationalize the off-by-one file loading
//and in the Portal, as of 8may07
//compiles and runs in Portal simulator
//begun delousing zatp.cpp for the Portal
//picks up 3 or the 5 files from the pwd=bin directory ... needs massive fixing!
//compiles and runs in the simulator for szg-1.0 now 
//but doesn't load the other molecules
//drawstars() has trouble, and the eyes,cube,cave,console structure
// main wants eyes==cube to compile in szg, lets try a hardcode to check
// eyes==console already compiles and works without the stars 
// I've commented out all uses of drawstars()
//continued Portal 7may07 
//gkf 30jul06, fixed a bunch of stuff, but the CAVE/CUBE version 
//still can't find its data. So, for now, we'll hardcode them here.
/***********************************************************************/
/*																	  **/
/*		computing frame on fly????	(would save memory requirements)  **/
/*		George's teapot code to rotate around newCenter 10/21/05********/
/*		Stuart's in-front-of-you code 10/14/05					********/
/********   Child of Entry for 2005 Scientific Visualization Challenge**/
/* 7/2/5 Need residueType[cAIndex] recorded in openpdbfile subroutine **/
/*	5/1/5deFa													********/
/*																	  **/
/*	install molecule scene.	vv										  **/
/*	hlight()				vv										  **/
/*	create some atomSpheres	vv										  **/
/*																	  **/
/*	install atomCenters[100][4] (4th coordinate for atom type, e.g., HYDRO)**/
/*																	  **/
/*			??attach center to sphere?? 							  **/
/*	could, by taking the north pole's "place"??						  **/
/*		?? need a whole row?										  **/
/*																	  **/
/*	10/26/05 initial run of loxodromic is stillborn	  **/
/*																	  **/
/*																	  **/
/*																	  **/
/* daughter of:														  **/
/*              Examples of Abelian Kleinian Groups:	    	*/
/*  Elementary Kleinian Groups are Finite Extensions of These	*/
/* Nov 22, 2004 modified to include Mac's			*/
/* June 14, 2004 modified for geowall		 		*/
/* June 8, 2004 Daughter of aKG.c, son of bishF04.c		*/
/****************************************************************/
/*		6/7/4 inserted loxodrome in the plane		*/
/*		6/10/4 disconnected ELLOX			*/
/****************************************************************/
/* Previous  work on akg.c:					*/
/* April 20, 1995 					 	*/
/* May 14, 1995							*/
/* Sept 25, 1995	Hyperbolic motion added		     	*/
/*  Extension of Mobius Transformation of H^2 into H^3 using quaternions*/
/* homotopy animation included March 30, 1995   */
/* Paul R. McCreary, Dept of Mathematics, Xula, all rights etc. */

/* (C) 2004 Board of Trustees Xavier Unversity of Louisiana */
/* e-mail pmccrear@xula.edu	*/

/* July 10th 2006 - David Bergman           */
/* Got atp.c to work in syzygy by reworking */
/* some lines of code - namely, including the precompiled header files,   */
/* fixing some typos that szygy was tripping up on and adding void to a number */
/* of functions that didn't actually return anything   */
/************************************************************************/
/*****************************************************************/
/****************************************************************/
/* Cross-platform definitions.                                  */
/****************************************************************/
#define LINUX 
//#define OSXMac I don't know how to handle these so I'm going to replace
// OSXMac with LINUX on Portal to see what happens
// #define WIN32

#define SZG

/**#define GeoWall**//*Obsolete*/
/*#define IRIX*/
/****************************************************************/
#include "arPrecompiled.h"  //must be the first line
#include "arMasterSlaveFramework.h"
// #include "arSZGClient.h"   zven.cpp doesn't have it in Portal


#include <stdlib.h>
#include <stdio.h>
#include <math.h>



#include <string.h>

#include "dawn.c"
#include <GL/glut.h>    //gkf 11jun07

#ifdef IRIX
#include <sys/time.h>           /* for the speedometer */
#include <glut.h>
#endif

#ifdef OSXMac
#include <sys/time.h>           /* for the speedometer */
#include <GLUT/glut.h>
#endif

#ifdef LINUX                  /* maybe this will work? */
#include <sys/time.h>
#endif

#ifdef WIN32
/* put this into the source code */
#pragma warning (disable:4305)  /* double-to-float  */
#pragma warning (disable:4101)	/* unreferenced local variable */
#pragma warning (disable:4056)	/* overflow in fp constant arithmetic */
#pragma warning (disable:4047)
#pragma warning (disable:4024)
#pragma warning (disable:4133)
#include <sys\timeb.h>          /* for the speedometer */
#include "arGlut.h"
//#include <GL\glut.h>
#define  M_PI  3.1415926535
#define  DATAPATH  g:/gfrancis/Data/zatp
#endif

#define  DATAPATH  /Users/georgefrancis/Data/molecules // this does not work nohow! 7may07portal 


#define  DTCONS     1           /* how fast autopilot is in console          */
#define  DTCAVE     10          /* how fast autopilot is in console          */
#define  BIGCAR     0.05        /* how BIG a coaster you want                */
#define  BIGCUBE    40.00       /* how BIG a cube you want                   */
#define Pi 3.1415926535

#define  ASDFG    1000        /* number of points computed                 */
#define  ZXCVB   (asdfg*2) /* number of development points              */

#define  LONGFLIGHT 500         /* how many contrails we record              */
#define  DGG 2*M_PI/(numbpts+1) /* Degrees?                                  */
#define  SLATEOFFSET 4.0        /* how far apart the development curves are  */

#define  MAX(x,y)       (((x)<(y))?(y):(x))
#define  MIN(x,y)       (((x)<(y))?(x):(y))
#define  ABS(u)         ((u)<0 ? -(u): (u))
#define  FOR(a,b,c)     for((a)=(b);(a)<(c);(a)++)
#define  DOT(p,q)       ((p)[0]*(q)[0]+(p)[1]*(q)[1]+(p)[2]*(q)[2])
#define  DOT4(p,q)      (DOT(p,q)+((p)[3]*(q)[3]))
#define  NRM(p)         sqrt(DOT((p),(p)))
#define  DG             M_PI/180
#define  S(u)           sin(u*DG)
#define  C(u)           cos(u*DG)
#define SCALE(A, B, C)  C[0]=A[0]*B; C[1]=A[1]*B; C[2]=A[2]*B;
#define MAG(A, B)       B = SQRT(DOT(A,A));
#define SQRT(u)         sqrt(u)
#define VSUB(A, B, C)   C[0]=A[0]-B[0]; C[1]=A[1]-B[1]; C[2]=A[2]-B[2];
#define VADD(A, B, C)   C[0]=A[0]+B[0]; C[1]=A[1]+B[1]; C[2]=A[2]+B[2];
#define CLAMP(x,u,v)    (x<u? u : (x>v ? v: x))
#define CROSS(A, B, C)  C[0]=A[1]*B[2]-A[2]*B[1];\
                        C[1]=A[2]*B[0]-A[0]*B[2];\
                        C[2]=A[0]*B[1]-A[1]*B[0];
#define FLYMODE         (1)
#define TURNMODE        (0)

/***artNscience: defining stuff*****************/

float circlePoints[1000][3];

/******* global variables ********/
		/* bio variables */

#define ALA 1
#define ARG 2
#define ASN 3
#define ASP 4
#define CYS 5
#define GLN 6
#define GLU 7
#define GLY 8
#define HIS 9
#define ILE 10
#define LEU 11
#define LYS 12
#define MET 13
#define PHE 14
#define PRO 15
#define SER 16
#define THR 17
#define TRP 18
#define TYR 19
#define VAL 20

       /* entry variables */
#define MARK 0 /* marker */
#define HYDR 1
#define CARB 2
#define NITR 3
#define OXYG 4
#define SULP 5
#define BKBN 9 /* backbone */
#define PHIL 10 /* hydrophilic */
#define PHOB 11 /* hydrophobic */
#define CHRG 12 /* charged side chains */
#define CHRGP 13 /* positively charged side chains */
#define CHRGN 14 /* negatively charged side chains */
#define ELSE 15 /* mics side chains */
#define NA2 16
#define CONINF 0 /* connectivity information */
#define RESTYP 1 /* residue type */
#define ATOTYP 2 /* atom type */
#define CHANID 3 /* chain id */
#define ATONUM 0 /* atom number */
#define RESNUM 1 /* residue number */
#define RESTYN 3 /* residue type (as number) */
#define NUMATOMS 20000 /* maximum number of atoms */
#define SPINENUM 30 /* num pnts btwn knots in cardinal spline */
#define CHAINUM 20 /* maximum num of chains */
#define LASTBBONE 5 /* number of backbone display options */
#define TAN 0
#define NOR 1
#define BIN 2
/*zxc*/
int testNum1,nnnn,openFileNum, bigMarker, begChainInx[100], showHydrophob=2;
int dummyNums[100][20],lastCent=NUMATOMS-1;
char dummyChars[/*NUMATOMS*/100][20][80], showMolName[80];
int nn1, nn2, nn3, residueType[1000], selectChains, whichChain=1, changeCent, axx;
int opaque, virg=1, flyForward, backBone, spineKnots, noConnections;
float atomRadius[10], rred, ggre, bblu, opaq, col[300][4], fastSpeed=4, slowSpeed=.5;
float SF,fac1=1,centerPntCorrect[3],tempaff[3], tempTan, tempNorm, tempBB;
float newCent[3],trans[3];
int targetAC, targetAC2,onetime,first1=0,fram1,onFly,chapmat;
int lastAtom, lastResidue, lastBB[CHAINUM], lastChain;
int chainCounter=1,atomCounter,jointCounter;
	/* zatp variables (shared flags)*/
int dfl2=1, changeScene=0, maxMol=4, changeMol=0;  //really 4 not 5!
       /* skel variables */
float gap, gap0=1.;                    /* kludge so that arguments() can set a default gap0 */
float lux[3]={1.,2.,3.};               /* world light direction vector                      */
float luxx[3];                         /* object space  direction vector                    */
float amb, pwr, lmb ;                       /* ambient fraction, pseudo-specular power           */
float mysiz,speed, torq, focal, wfar;  /* console navigation variables                      */
int win= 2;                           /* 2 full screen, use 0 for demand sized             */
unsigned int BUT, XX, YY,SHIF;         /* used in chaptrack gluttery                        */
int xt,yt;                             /* once was xt,yt,xm,ym for viewportery              */
int mode,morph,msg,binoc;              /* pretty global                                     */
int th0, th1, dth, ta0, ta1, dta;      /* torus parameters                                  */

float t=0, temp2, temp3;               /* we hope                                           */
int ii, jj, kk;  float tmp, temp;      /* saves gray hairs later                            */
float aff[16], aff2[16], starmat[16], mat[16];   /* places for the world                              */
float mat1[4][4],mat2[4][4],mat3[4][4];
float nose;                            /* to eye distance in console                        */
int hasnumber, number, decimal, sign ; /* for Stuart's number parser                        */
float dot1, dot2, dot3;
float tempvec[3];

int showcube;

       /* bishop variables */
int mousePosition[3];                  /* the mouse position for picking          */
int friz;                              /* friz is back                            */
int whiteBG;			/* white or black background */
#define X 0	/* x-coordinate or auxiliary "type" of frame or still */
#define Y 1
#define Z 2
#define W 3
#define Pi 3.1415926535
#define NDP 20	/* NDP = number of points in "dot" */
#define NPP 4*2	/* NPP = 2*number of points across path */
#define  GETS(pt1,pt2)  FOR(inx,0,3){pt1[inx]=pt2[inx];}
#define  GETS2D(pt1,pt2)  FOR(inx,0,2){pt1[inx]=pt2[inx];}
int inx, III,scene;
float dddd,aaa=.06,bbb=1,hhhh=-.2;
float Molecule[100][20][NDP][NDP][3];
float HexVerts[10][30][3],PentVerts[10][30][3];
float AdenineVerts[100][30][3], ThymineVerts[100][30][3],GuanineVerts[100][30][3],CytosineVerts[100][30][3];

/** asdf **/
float speedMed;
int XXX=0, YYY=1, ZZZ=2;	/* coordinates */
FILE* openptr1;
/**akgeowall**/
#define LOXODROMIC 0
#define ELLIPTIC 1
#define KLEIN4 2
/*#define ELLOX 3 deactivated*/
#define PARABOLIC 3
#define MOLECULE  4
#define SPARCE 1
#define VERBOSE 2
#define loxPlaneColor  glColor3f(0.2,1.0,1.0)
#define loxAxisColor  glColor3f(1.0,1.0,0.2)
#define luneColor  glColor3f(1.0,0.2,1.0)
#define NAP 57	/* NAP = number of points along central axis */
#define NA 51
#define NR 4  /* number of sets of rings */
#define NL 10  /* number of lunes */
#define TH  5  /* number of stablized layers to lune */
#define NHAP 200   /* number of points on each circle of pencil */
#define NH  8   /* number of circles in the hyperbolic pencil */
#define NLATS 10  /* number of latitudes */
#define NLONGS 10  /* number of longitudes */
#define AAAA 0
#define BBBB 1
#define CCCC 2
#define DDDD 3
float test1[3];
int geo=0,numIterates,cMax,horo,sceneNum,loxCount,ellCount,thisOne,slideInx,lMax=10;
int pla,justStarted,mov,rot,klein2,numLayers=1,pauseCounter,pauseMax;
float one[3]={1,0,0},negOne[3]={-1,0,0},origin[3]={0,0,0},cf[100][4][3],radd1,radd2;
float paraPlane[200][200][3],fPLats[NLATS][NAP][3],fPLongs[NLONGS][NAP][3],
	frameHypPencil[NH][NHAP][3],hypPencil[NH][NHAP][3],sLats[NLATS][NAP][3],
	sLongs[NLONGS][NAP][3],pLats[NLATS][NAP][3],pLongs[NLONGS][NAP][3],
	tetra[2][4][3][3][NAP][3],frameTetra[2][4][3][3][NAP][3],baseVertices[3][3],
	slidePnt[10][3];
/*end akgeowall*/
/*frames*/
#define NUMCPnts 1000
float development1[2][NUMCPnts], development2[2][NUMCPnts];         /* developments of the frames */
float curvefrenetframe[3][NUMCPnts],curvebishopframe[3][NUMCPnts];
float frenetcurve[3][NUMCPnts],bishopcurve[3][NUMCPnts];
int whichFrame;
       /* bishop variables */
int numbpts, paramlength;                           /* variable number of points on the curves */
float initpos[3],initframe[9];
/*end of frames*/

/**********************CAVEISH****************************************/
#ifndef OSXMac
#include <malloc.h>     /* for the shared vars */
#endif

#ifdef CAVE
#include <cave_ogl.h>
#include <cave.macros.h>
#endif  /* CAVE */

int caveyes ;  /* read cave...yes! */

/*Shared memory variables as per Stuart Levy 1994 */
/*********************************************************************/
/* This shared variable structure is for the CAVE.  I truly don't    */
/*    know if I should be sharing everything that I do.  As a general*/
/*    rule, if I need a global variable I put it here. ~baf 31jul02  */
/*********************************************************************/
struct share_var {
	/*akgeowall*/
   int s_bbIndices[CHAINUM][NUMATOMS/3];
   float s_controlPnts[NUMATOMS][3], s_backbone[CHAINUM][NUMATOMS*SPINENUM][3];
   char s_atomInfoStr[NUMATOMS][4][80];
   int s_atomInfoNum[NUMATOMS][4], s_alphaCIndices[NUMATOMS];
   float s_drawSph[NUMATOMS][NA2+1][NA2+1][3], s_aCent[NUMATOMS][4];
   float s_aSphere[NA2+1][NA2+1+1][3];
   int s_loxLength,s_loxLengthPerm;
   float s_loxCurvePoints[2][1000][3], s_frameLoxCurvePoints[2][1000][3],s_frameLune[NL][TH][NA][NAP][3],s_lune[NL][TH][NA][NAP][3];
   float s_frameRing[NL][NR][10][NAP][3],s_ring[NL][20][10][5][NAP][3],s_frameAxes[NL][NAP][3],s_axes[NL][NAP][3];
   float s_equator[NL][16][NAP][3];
   float s_twist,s_loxTwist,s_multiplierK,s_fixedPnts[10][2][3],s_theta[TH],s_infToF1[6][4][3],s_f1ToInf[6][4][3];
   int s_band,s_loxPath,s_lastLune,s_order,s_sceneLimit,s_movieCounter,s_sceneType;
   /*end akgeowall*/
   float  s_siz;            /* final scaling factor before projection     */
   int s_wnd;               /* wand or maus flag                          */
   int s_gnd;               /* background color                           */
   float s_luxx[3];         /* current light direction                    */
   int s_mauspaw ;          /* binary for mouse button state              */
   int s_th0 ;              /* since the autotymer changes these          */
   int s_th1 ;              /* they need to be shared                     */
   int s_ta0 ;
   int s_ta1 ;
   float s_gap ;
   float  s_aa ;            /* Karan .... document these                  */
   float  s_bb ;
   float  s_cc ;
   float s_sizz;
   int  s_dtt ;
   int  s_tt ;
   GLfloat s_aff[16], s_starmat[16], s_aff2[16]; /* affine matrices for object, stars */
   GLfloat s_const_aff[16], s_const_aff2[16];          /* constant affine matrix            */


   float s_contrail[LONGFLIGHT][16];

   int s_showcontrail;     /* flag to show contrail                */
   int s_showslate;        /* flag to show slate                   */
   int s_cockpit;          /* flag to view from cockpit            */
   int s_showbishop;       /* flag to show the frenet/bishop frame */
   int s_autopilot;        /* flag to turn on autopilot            */

} *s_var;
	/*akgeowall*/
#define bbIndices (s_var->s_bbIndices)
#define controlPnts (s_var->s_controlPnts)
#define backbone (s_var->s_backbone)
#define atomInfoStr (s_var->s_atomInfoStr)
#define alphaCIndices (s_var->s_alphaCIndices)
#define atomInfoNum (s_var->s_atomInfoNum)
#define aCent  (s_var->s_aCent)
#define aSphere (s_var->s_aSphere)
#define loxLength (s_var->s_loxLength)
#define loxLengthPerm	(s_var->s_loxLengthPerm)
#define loxCurvePoints (s_var->s_loxCurvePoints)
#define frameLoxCurvePoints (s_var->s_frameLoxCurvePoints)
#define frameLune	(s_var->s_frameLune)
#define lune		(s_var->s_lune)
#define frameRing	(s_var->s_frameRing)
#define ring		(s_var->s_ring)
#define frameAxes	(s_var->s_frameAxes)
#define axes		(s_var->s_axes)
#define equator		(s_var->s_equator)
#define twist		(s_var->s_twist)
#define loxTwist	(s_var->s_loxTwist)
#define multiplierK	(s_var->s_multiplierK)
#define fixedPnts	(s_var->s_fixedPnts)
#define theta		(s_var->s_theta)
#define infToF1		(s_var->s_infToF1)
#define f1ToInf		(s_var->s_f1ToInf)
#define band		(s_var->s_band)
#define loxPath		(s_var->s_loxPath)
#define lastLune	(s_var->s_lastLune)
#define order		(s_var->s_order)
#define sceneLimit	(s_var->s_sceneLimit)
#define movieCounter (s_var->s_movieCounter)
#define sceneType	(s_var->s_sceneType)
	/*end akgeowall*/
#define siz       (s_var->s_siz)
#define wnd       (s_var->s_wnd)
#define gnd       (s_var->s_gnd)  /* use only if printframe curve is present */
#define luxx      (s_var->s_luxx)
#define mauspaw   (s_var->s_mauspaw)
#define sizz     (s_var->s_sizz)
#define th0       (s_var->s_th0)
#define th1       (s_var->s_th1)
#define ta0       (s_var->s_ta0)
#define ta1       (s_var->s_ta1)
#define gap       (s_var->s_gap)
#define aa        (s_var->s_aa)
#define bb        (s_var->s_bb)
#define cc        (s_var->s_cc)
#define dtt       (s_var->s_dtt)
#define tt        (s_var->s_tt)
#define aff       (s_var->s_aff)
#define aff2       (s_var->s_aff2)
#define starmat   (s_var->s_starmat)
#define const_aff (s_var->s_const_aff)
#define const_aff2 (s_var->s_const_aff2)

#define contrail (s_var->s_contrail)
#define showslate      (s_var->s_showslate)
#define cockpit        (s_var->s_cockpit)
#define showbishop     (s_var->s_showbishop)
#define autopilot      (s_var->s_autopilot)

/**********************************************************************/
/* function prototypes - fixes functions being defined out of order   */
/*    so that I can group functions according to their purpose.  This */
/*    is only a partial list.                                         */
/**********************************************************************/
void drawall(void);
void drawall2(void);
void drawallllll(void);
void openpdbfile(void);
void computeCardinalSpline2(float P0,float P1,float P2,float P3,float t);
/**********************************************************************/
/* getmem - allocates the memory initially.                           */
/**********************************************************************/
void getmem(void){
    /* remember, you need (struct share_var*) for C++ compiler */
    if(caveyes)
#ifdef CAVE
    s_var = (struct share_var *) CAVEMalloc(sizeof(struct share_var))
    /* no ; here ! */
#endif /* CAVE */
     ; /* but ;  here */
   else /* CONSOLE */
     s_var =(struct share_var *) malloc(sizeof(struct share_var));
    if(s_var==NULL){fprintf(stderr,"No room to share! %x \n",s_var); exit(1);}

}
/******************end*CAVEISH*****************************************/

/******************gesture**********************************/
int gesture =0; 
float eps =.002;   //gain factor for rotations
float del =.02;   //gain factor for translations
int resethand =0;  
/******************gesture**********************************/


/**********************************************************************/
/* autotymer - illiView animator                                      */
/**********************************************************************/
void autotymer(int rreset){ /* cheap animations */
#ifdef IRIX
#define  TYME(cnt,max,act) {static cnt; if(first)cnt=max; else\
                            if(cnt?cnt--:0){ act ; goto Break;}}
  static first = 1;  /* the first time autymer is called */
  if(rreset)first=1;  /* or if it is reset to start over  */
  TYME( shrink , 150,th0++;th1--;ta0++;ta1--)
  TYME( pause  , 20,          )
  TYME( grow   , 150,th0--;th1++;ta0--;ta1++)
  TYME( dwell  , 30,            )
  TYME(finish  , 1 , first = 1  )
  first = 0;
  Break:   ;   /* yes Virginia, C has gotos */
#endif
}

void recordAFF(void){ /** print out affine matrix with a keystroke:
                used to record aff matrix of any set position **/
printf("aff[0*4+0]=%g;",aff[0*4+0]);
printf("aff[0][1]=%g;",aff[0*4+1]);
printf("aff[0][2]=%g;",aff[0*4+2]);
printf("aff[0][3]=%g;\n",aff[0*4+3]);
printf("aff[1][0]=%g;",aff[1*4+0]);
printf("aff[1][1]=%g;",aff[1*4+1]);
printf("aff[1][2]=%g;",aff[1*4+2]);
printf("aff[1][3]=%g;\n",aff[1*4+3]);
printf("aff[2][0]=%g;",aff[2*4+0]);
printf("aff[2][1]=%g;",aff[2*4+1]);
printf("aff[2][2]=%g;",aff[2*4+2]);
printf("aff[2][3]=%g;\n",aff[2*4+3]);
printf("aff[3][0]=%g;",aff[3*4+0]);
printf("aff[3][1]=%g;",aff[3*4+1]);
printf("aff[3][2]=%g;",aff[3*4+2]);
printf("aff[3][3]=%g;\n\n",aff[3*4+3]);
}

void recordAFF2(void){ /** print out affine matrix with a keystroke:
                used to record aff matrix of any set position **/
printf("aff2[0*4+0]=%g;",aff2[0*4+0]);
printf("aff2[0][1]=%g;",aff2[0*4+1]);
printf("aff2[0][2]=%g;",aff2[0*4+2]);
printf("aff2[0][3]=%g;\n",aff2[0*4+3]);
printf("aff2[1][0]=%g;",aff2[1*4+0]);
printf("aff2[1][1]=%g;",aff2[1*4+1]);
printf("aff2[1][2]=%g;",aff2[1*4+2]);
printf("aff2[1][3]=%g;\n",aff2[1*4+3]);
printf("aff2[2][0]=%g;",aff2[2*4+0]);
printf("aff2[2][1]=%g;",aff2[2*4+1]);
printf("aff2[2][2]=%g;",aff2[2*4+2]);
printf("aff2[2][3]=%g;\n",aff2[2*4+3]);
printf("aff2[3][0]=%g;",aff2[3*4+0]);
printf("aff2[3][1]=%g;",aff2[3*4+1]);
printf("aff2[3][2]=%g;",aff2[3*4+2]);
printf("aff2[3][3]=%g;\n\n",aff2[3*4+3]);
}

void buildFromAFF(float R[4][4]){
	int i,j;
	FOR(i,0,4){
		FOR(j,0,4){
			R[i][j]=aff[i*4+j];
		}
	}
}
void assignAFF2(float R[4][4]){
	int i,j;
	FOR(i,0,4){
		FOR(j,0,4){
			aff2[i*4+j]=R[i][j];
		}
	}
}
void multMatrices(float R[4][4],float M[4][4],float N[4][4]){
	int i,j;
	FOR(i,0,4){
		FOR(j,0,4){
			R[i][j]=M[i][0]*N[0][j]+M[i][1]*N[1][j]+M[i][2]*N[2][j]+M[i][3]*N[3][j];
		}
	}
}
void invertMatrix(float R[4][4], float M[4][4]){
	int i,j;
	FOR(i,0,3){
		FOR(j,0,3){
			R[i][j]=M[j][i];
		}
		R[3][i]=-M[3][i];
	}R[3][3]=M[3][3];
}
void checkProcess(void){
	buildFromAFF(mat1);
	invertMatrix(mat2,mat1);
	/**multMatrices(mat3,mat1,mat2);
	assignAFF2(mat3);/**/

/**/assignAFF2(mat2);/**/

	recordAFF2();

}
void pushit(float *a){
  a[0] = aff[0]*newCent[X]+aff[4]*newCent[Y]+aff[8]*newCent[Z]/**/+aff[12]/**/;
  a[1] = aff[1]*newCent[X]+aff[5]*newCent[Y]+aff[9]*newCent[Z]/**/+aff[13]/**/;
  a[2] = aff[2]*newCent[X]+aff[6]*newCent[Y]+aff[10]*newCent[Z]/**/+aff[14]/**/;
}
void pushit2(float *a){
  a[0] = aff2[0]*newCent[X]+aff2[4]*newCent[Y]+aff2[8]*newCent[Z]/**/+aff2[12]/**/;
  a[1] = aff2[1]*newCent[X]+aff2[5]*newCent[Y]+aff2[9]*newCent[Z]/**/+aff2[13]/**/;
  a[2] = aff2[2]*newCent[X]+aff2[6]*newCent[Y]+aff2[10]*newCent[Z]/**/+aff2[14]/**/;
}
void pushit3(float *a){
  a[0] = aff2[0]*newCent[X]+aff2[4]*newCent[Y]+aff2[8]*newCent[Z]+aff2[12]+aff[12];
  a[1] = aff2[1]*newCent[X]+aff2[5]*newCent[Y]+aff2[9]*newCent[Z]+aff2[13]+aff[13];
  a[2] = aff2[2]*newCent[X]+aff2[6]*newCent[Y]+aff2[10]*newCent[Z]+aff2[14]+aff[14];
}
/*******************
#define skW .05

float ax[21][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
        {1.3,0,0},{1.1,0,skW},{1.1,0,-skW},{1.1,skW,0},{1.1,-skW,0},
        {0,1.3,0},{skW,1.1,0},{-skW,1.1,0},{0,1.1,skW},{0,1.1,-skW},
        {0,0,1.3},{skW,0,1.1},{-skW,0,1.1},{0,skW,1.1},{0,-skW,1.1}};
/*
*       SUBROUTINE "axes()" GOES IN DRAWIT ROUTINE
*/
/*****************
drawAxes(){
        bgnline();cpack(0xff0000ff);v3f(ax[0]);v3f(ax[1]);endline();
        bgnline();cpack(0xff00ff00);v3f(ax[2]);v3f(ax[3]);endline();
        bgnline();cpack(0xffff0000);v3f(ax[4]);v3f(ax[5]);endline();
        bgntmesh();cpack(0xff0000ff);v3f(ax[6]);v3f(ax[7]);swaptmesh();
        v3f(ax[9]);swaptmesh();v3f(ax[8]);swaptmesh();v3f(ax[10]);swaptmesh();v3f(ax[7]);endtmesh();
        bgntmesh();cpack(0xff00ff00);v3f(ax[11]);v3f(ax[12]);swaptmesh();v3f(ax[14]);swaptmesh();v3f(ax
[13]);swaptmesh();v3f(ax[15]);swaptmesh();v3f(ax[12]);endtmesh();
        bgntmesh();cpack(0xffff0000);v3f(ax[16]);v3f(ax[17]);swaptmesh();v3f(ax[19]);swaptmesh();v3f(ax
[18]);swaptmesh();v3f(ax[20]);swaptmesh();v3f(ax[17]);endtmesh();
}

/**********************************************************************/
/**********************************************************************/

/****************Stuart's findscreen*********************/

/* These values are constants during the time a given frame is being drawn.
 * Update them at beginning of each frame by calling findscreen().
 */
float pixelsperradian;
float forward[3];
float eyepoint[3];


/* dst[] and src[] are 3-vectors; T is a 4x4 matrix */
void vtfmvector2( float dst[3], float src[3], float T[4*4] )
{
  int i;
  for(i = 0; i < 3; i++)
    dst[i] = src[0]*T[0*4+i] + src[1]*T[1*4+i] + src[2]*T[2*4+i];
}

/* Invert a matrix, assuming it's a Euclidean isometry
 * plus possibly uniform scaling.
 */
void eucinv( float Tdst[4*4], float Tsrc[4*4] )
{
  int i, j;
  float s = DOT(Tsrc, Tsrc);	/* measure scaling by looking at first row */
  float trans[3];

  for(i = 0; i < 3; i++) {
    for(j = 0; j < 3; j++)
	Tdst[i*4+j] = Tsrc[j*4+i] / s;
    Tdst[i*4+3] = 0;
  }
  vtfmvector2( trans, &Tsrc[4*3+0], Tdst );
  Tdst[3*4+0] = -trans[0];
  Tdst[3*4+1] = -trans[1];
  Tdst[3*4+2] = -trans[2];
  Tdst[3*4+3] = 1;
}

/* Reads OpenGL matrices; determines...
 *   eyepoint (in world coordinates)
 *   forward  (unit vector in world coords, perp. to screen)
 *   pixelsperradian (scale factor: angular size to size on screen)
 */
void findscreen()
{
  float Tw2c[4*4];	// world-to-eye (MODELVIEW) matrix
  float Tproj[4*4];	// eye-to-normalized-device (PROJECTION) matrix
  int xywh[4],i;		// screen (well, viewport) size

  float Tc2w[4*4];

  glGetFloatv( GL_MODELVIEW_MATRIX, Tw2c );
  eucinv( Tc2w, Tw2c );


  glGetFloatv( GL_PROJECTION_MATRIX, Tproj );
  glGetIntegerv( GL_VIEWPORT, xywh );
  pixelsperradian = (.5*xywh[2] * Tproj[0*4+0]);

  /* Construct a "forward" vector in object coords,
   * for measuring distance from camera plane.
   * Note camera looks toward its -Z axis not +Z!
   */
  forward[0] = -Tc2w[2*4+0]*fac1;
  forward[1] = -Tc2w[2*4+1]*fac1;
  forward[2] = -Tc2w[2*4+2]*fac1;

  /* and extract the world-coordinate eyepoint from Tc2w */
  eyepoint[0] = Tc2w[3*4+0];
  eyepoint[1] = Tc2w[3*4+1];
  eyepoint[2] = Tc2w[3*4+2];

  /**/if(changeCent){
	  	FOR(i,0,3){
			newCent[i]=aCent[lastCent][i]=forward[i]+eyepoint[i];
	  	}
	  }/**/

}

/* Compute the screen-radius in pixels of
 * a sphere at position worldpos, of given 3-D radius.
 * Assumes that findscreen() has been called for this frame already.
 */
float pixelsize( float *worldpos, float radius )
{
  float there[3];	/* vector from worldpos to eyepoint */
  float distance;

  VSUB( worldpos,eyepoint, there );
  distance = NRM( there );
  return pixelsperradian * radius / distance;
}
/*******************End of Stuart's findscreen*************************/


/*****************************************************************/
/*******  COMPLEX ARITHMETIC  ************************************/
/*****************************************************************/

void add(float rp[3], float p1[3], float p2[3]){
	rp[X]=p1[X]+p2[X];
	rp[Y]=p1[Y]+p2[Y];
}
void sub(float rp[3], float p1[3], float p2[3]){
	rp[X]=p1[X]-p2[X];
	rp[Y]=p1[Y]-p2[Y];
}
void mul(float rp[3], float p1[3], float p2[3]){
	rp[X]=p1[X]*p2[X]-p1[Y]*p2[Y];
	rp[Y]=p1[X]*p2[Y]+p1[Y]*p2[X];
}
void div2(float rp[3], float p1[3], float p2[3]){
	rp[X]=(p1[X]*p2[X]+p1[Y]*p2[Y])/(p2[X]*p2[X]+p2[Y]*p2[Y]);
	rp[Y]=(-p1[X]*p2[Y]+p1[Y]*p2[X])/(p2[X]*p2[X]+p2[Y]*p2[Y]);
}
void cnj(float rp[3], float p1[3])/*float rp[3], p1[3];*/{
	rp[X]=p1[X];
	rp[Y]=-p1[Y];
}
float absValue2DD(float pt[3]){
        return (pow(pow(pt[X],2.f)+pow(pt[Y],2.f),.5f));
}
void sqr(float rp[3], float p1[3]){
	float r, thet;
	r=pow(absValue2DD(p1),.5f);
        if(p1[Y]>0) thet=acos(p1[X]/absValue2DD(p1));
        else thet=2*Pi-acos(p1[X]/absValue2DD(p1));
	rp[X]=r*cos(thet/2);
	rp[Y]=r*sin(thet/2);
}
void slide(float resPnt[3], float pnt[3], float transPnt[3]){
	int i;
	FOR(i,0,3){
		resPnt[i]=pnt[i] + transPnt[i];
	}
}
void multMatrix(float res[4][3],float m1[4][3],float m2[4][3]){
	int i;
	float temp[4][3];
	FOR(i,0,4){GETS(temp[i],origin);}
	mul(temp[0],m1[AAAA],m2[AAAA]);
	mul(temp[1],m1[BBBB],m2[CCCC]);
	add(res[AAAA],temp[0],temp[1]);
	mul(temp[0],m1[AAAA],m2[BBBB]);
	mul(temp[1],m1[BBBB],m2[DDDD]);
	add(res[BBBB],temp[0],temp[1]);
	mul(temp[0],m1[CCCC],m2[AAAA]);
	mul(temp[1],m1[DDDD],m2[CCCC]);
	add(res[CCCC],temp[0],temp[1]);
	mul(temp[0],m1[CCCC],m2[BBBB]);
	mul(temp[1],m1[DDDD],m2[DDDD]);
	add(res[DDDD],temp[0],temp[1]);
}
/**spherical geometry *****/
void rotateT(float resultPt[3], float point[3], float angle){
        float initialAngle, r;
        r=absValue2DD(point);
        if(r!=0){
                if(point[Y]>0) initialAngle=acos(point[X]/r);
                else initialAngle=2*Pi-acos(point[X]/r);
                resultPt[X]=r*cos(angle+initialAngle);
                resultPt[Y]=r*sin(angle+initialAngle);
                resultPt[Z]=point[Z];
        }
        else 	GETS(resultPt,point);
}
/* pE(resultPoint,(x,y,t),(a,b,c,d)): Poincare extension of the mobius transformation,
*****             (az + b)/(cz + d), where a, b, c, & d are all complex
*****Summary: a point {x,y,t} goes to {Re[pE(x,y,t)],Im[pE(x,y,t)],J[pE(x,y,t)]}*/
/***cf[10][4][3]*/
void pE(float rp[3], float pnt[3], float c[4][3]){
	float zz[3]={0,0,0},t,
			t1[3]={1,1,1},t2[3]={1,1,1},t3[3]={1,1,1},t4[3]={1,1,1},
			denominator;
	GETS2D(zz,pnt);t=pnt[Z];
	mul(t1,c[CCCC],zz);/* c*z */
	add(t2,t1,c[DDDD]);/* c*z + d */
	cnj(t4,t2); /* conj(c*z + d) SAVE for numerator*/
	mul(t1,t2,t4);/* |c*z + d|^2 SAVE for denom */

	cnj(t2,c[CCCC]);/* conj(c) */
	mul(t3,c[CCCC],t2); /* |c|^2 */

	denominator = t1[X] + t3[X]*t*t;
		/* END COMPUTATION SECTION, but save t4 */

	mul(t1,c[AAAA],c[DDDD]);
	mul(t2,c[BBBB],c[CCCC]);
	sub(t3,t1,t2);/* (ad-bc) */
	cnj(t1,t3);
	mul(t2,t3,t1);/* |ad-bc|^2 */

	rp[Z] = pow(t2[X],.5f)*t/denominator;
		/* END COMPUTATION (j COMPONENT)*/

	/********** complex component  ************************/
	mul(t1,c[AAAA],zz);/* az */
	add(t2,t1,c[BBBB]);/* az+b */
	mul(t3,t2,t4);/*t4 from above->(az+b)*conj(cz+d)*/
	cnj(t1,c[CCCC]);
	mul(t2,c[AAAA],t1);/* a*conj(c) */
	/* t2=a*conj(c), t3=az+b)*conj(cz+d) */
	rp[X] = (t3[X] + t2[X]*t*t)/denominator;
	rp[Y] = (t3[Y] + t2[Y]*t*t)/denominator;

}
void stereoProjection(float res[3],float point[3]){
	res[X]= 2*point[X]/(pow(point[X],2)+pow(point[Y],2)+1);
        res[Y]= 2*point[Y]/(pow(point[X],2)+pow(point[Y],2)+1);
        res[Z]=(pow(point[X],2)+pow(point[Y],2)-1)/(pow(point[X],2)+pow(point[Y],2)+1);
}
void inverseStereoProjection(float res[3],float point[3]){
        res[X] = .5*(point[X]*( (point[Z] + 1)/(1 - point[Z]) + 1) );
        res[Y] = .5*(point[Y]*( (point[Z] + 1)/(1 - point[Z]) + 1) );
}
void newCoefficients(int whichOne,float fix1[3], float fix2[3], float multiplier[3]){
	int n;
	float mu[3], imu[3], t1[3], t2[3], t3[3];
	sqr(mu,multiplier); /* mu = SQRT(K) */
	div2(imu,one,mu);
	n=whichOne;
	mul(t1,fix1,mu);
	mul(t2,fix2,imu);
	sub(cf[n][AAAA],t1,t2);
	sub(t1,imu,mu);
	mul(t2,fix1,t1);
	mul(cf[n][BBBB],fix2,t2);
	sub(cf[n][CCCC],mu,imu);
	mul(t1,fix1,imu);
	mul(t2,fix2,mu);
	sub(cf[n][DDDD],t1,t2);
}
void newFaces(void){/*tetra[t][f][v][a][p]*/
	int i, j, k, v/*ertices*/, a/*rcs*/, p/*oints*/, f/*ace*/;
	float  two[3]={2,0,0}, multiplier[2], te[3],tm[4][3],tmp2[4][3],
		slideP[3];
	FOR(i,0,4){GETS(tm[i],origin);GETS(tmp2[i],origin);}
	slideP[X]=radd1;slideP[Y]=radd1*pow(3.f,.5f);slideP[Z]=0;
	multiplier[X]=cos(-Pi/3);multiplier[Y]=sin(-Pi/3);
	newCoefficients(4,baseVertices[0],baseVertices[1],multiplier);
	FOR(i,0,3){slidePnt[0][i]=baseVertices[0][i]-baseVertices[2][i];}
	FOR(i,0,3){slidePnt[1][i]=baseVertices[2][i]-baseVertices[0][i];}
	FOR(i,0,3){slidePnt[2][i]=baseVertices[1][i]-baseVertices[2][i];}
	FOR(i,0,3){slidePnt[3][i]=baseVertices[2][i]-baseVertices[1][i];}
	FOR(i,0,3){
		FOR(j,0,3){
			FOR(k,0,NAP){
				pE(tetra[0][1][i][j][k],
					tetra[0][0][i][j][k],cf[4]);
			}
		}
	}
	FOR(i,0,3){
		FOR(j,0,3){
			FOR(k,0,NAP){
				rotateT(te,tetra[0][1][i][j][k],2*Pi/3);
				slide(tetra[0][2][i][j][k],te,one);
				rotateT(te,tetra[0][1][i][j][k],4*Pi/3);
				slide(tetra[0][3][i][j][k],te,slideP);
			}
		}
	}
	FOR(f,0,4){
		FOR(v,0,3){
			FOR(a,0,3){
				FOR(p,0,NAP){
		  pE(tetra[1][f][v][a][p],tetra[0][f][v][a][p],cf[4]);
				}
			}
		}
	}

}
void drawAxes(){
	float orig[3],bend[3], fac2=10;

       orig[X]=0;orig[Y]=0;orig[Z]=0;
       glLineWidth(5);
	glColor3f(1.0,0.0,0.0);bend[X]=1*fac2;bend[Y]=0;bend[Z]=0;
       glBegin(GL_LINE_STRIP);
		glVertex3fv(orig);glVertex3fv(bend);
       glEnd();
	glColor3f(0.0,1.0,0.0);bend[X]=0;bend[Y]=1*fac2;bend[Z]=0;
       glBegin(GL_LINE_STRIP);
		glVertex3fv(orig);glVertex3fv(bend);
       glEnd();
	glColor3f(0.0,0.0,1.0);bend[X]=0;bend[Y]=0;bend[Z]=1*fac2;
       glBegin(GL_LINE_STRIP);
		glVertex3fv(orig);glVertex3fv(bend);
       glEnd();
}
/**artNscience: construction of objects***********/
void drawTheCircle(void){
	int ii, jj, numPnts;
	numPnts=51;
      glLineWidth(5);
	glColor3f(1.0,1.0,1.0);
       glBegin(GL_LINE_STRIP);
	FOR(ii,0,numPnts){
		glVertex3fv(circlePoints[ii]);
	}
       glEnd();
}
void constructACircle(void){
	int ii, jj, numPnts;
	float centerOfCircle[3], circleRadius;
	circleRadius=2;
	centerOfCircle[X]=1;
	centerOfCircle[Y]=.5;
	centerOfCircle[Z]=0;
	numPnts=51;
	FOR(ii,0,numPnts){
		circlePoints[ii][X]=cos(ii*2*Pi/numPnts);
		circlePoints[ii][Y]=sin(ii*2*Pi/numPnts);
		circlePoints[ii][Z]=0;
	}
}/*****************/
void constructTetrahedron(void){
	int i, j, k, p;
	float theta2;
	float epsilon;
	float basicEdge[NAP][3];
	float te[3];
	float slideTemp[4][3];
	radd1=.5;
	epsilon=.0001;
	FOR(i,0,NAP){
		theta2=(Pi/2-epsilon)*(float)i/(NAP-1) + epsilon;
		basicEdge[i][Y]=0;
		basicEdge[i][X]=radd1*cos(theta2);
		basicEdge[i][Z]=radd1*sin(theta2);
	}
	slideTemp[0][X]=radd1;slideTemp[0][Y]=0;slideTemp[0][Z]=0;
	slideTemp[1][X]=radd1/2;slideTemp[1][Y]=radd1*pow(3.f,.5f)/2;slideTemp[1][Z]=0;
	slideTemp[2][X]=3*radd1/2;slideTemp[2][Y]=radd1*pow(3.f,.5f)/2;slideTemp[2][Z]=0;
	FOR(i,0,NAP){
		rotateT(te,basicEdge[i],Pi);
		slide(tetra[0][0][0][0][i],te,slideTemp[0]);
		rotateT(te,basicEdge[i],Pi+Pi/3);
		slide(tetra[0][0][0][2][i],te,slideTemp[1]);
		rotateT(te,basicEdge[i],5*Pi/3);
		slide(tetra[0][0][1][0][i],te,slideTemp[2]);
		slide(tetra[0][0][1][2][i],basicEdge[i],slideTemp[0]);
		rotateT(te,basicEdge[i],Pi/3);
		slide(tetra[0][0][2][0][i],te,slideTemp[1]);
		rotateT(te,basicEdge[i],2*Pi/3);
		slide(tetra[0][0][2][2][i],te,slideTemp[2]);
	}
	radd2=2./3.*radd1*pow(3.f,.5f);
	FOR(i,0,NAP){
		theta2=(Pi/2-epsilon)*(float)i/(NAP-1) + epsilon;
		basicEdge[i][Y]=0;
		basicEdge[i][X]=radd2*cos(theta2);
		basicEdge[i][Z]=radd2*sin(theta2);
	}
	slideTemp[0][X]=radd2*pow(3.f,.5f)/2;slideTemp[0][Y]=radd2/2;slideTemp[0][Z]=0;
	slideTemp[1][X]=2*radd1;slideTemp[1][Y]=0;slideTemp[1][Z]=0;
	slideTemp[2][X]=radd1;slideTemp[2][Y]=radd1*pow(3.f,.5f);slideTemp[2][Z]=0;
	FOR(i,0,NAP){
		rotateT(te,basicEdge[i],Pi+Pi/6);
		slide(tetra[0][0][0][1][i],te,slideTemp[0]);
		rotateT(te,tetra[0][0][0][1][i],2*Pi/3);
		slide(tetra[0][0][1][1][i],te,slideTemp[1]);
		rotateT(te,tetra[0][0][0][1][i],4*Pi/3);
		slide(tetra[0][0][2][1][i],te,slideTemp[2]);
	}
	newFaces();
}
void grids(){  /*** All lunes ***************/
	int i, j, k, l, m, n, nn, torder, p, r;
	float a[3], b[3],c[3],d[3],f1[3],f2[3],denom,fac,begFac,multiplier,begFacEL,temp[3];
	multiplier=1/.9;
	begFac=pow(1.f/multiplier,NA/2.f-1.f);
	FOR(nn,0,2){ /** define mobius transformations to move 0 & INF to f2 & f1 ****/
		GETS(f1,fixedPnts[nn][0]); GETS(f2,fixedPnts[nn][1]);
		denom=pow(pow(-f1[X]+f2[X],2.f)+pow(-f1[Y]+f2[Y],2.f),.25f);
	  	infToF1[nn][AAAA][X]=-f1[X]/denom;infToF1[nn][BBBB][X]=f2[X]/denom;infToF1[nn][CCCC][X]=-1/denom;infToF1[nn][DDDD][X]=1/denom;
	  	infToF1[nn][AAAA][Y]=-f1[Y]/denom;infToF1[nn][BBBB][Y]=f2[Y]/denom;infToF1[nn][CCCC][Y]=0;	 infToF1[nn][DDDD][Y]=0;
	  	f1ToInf[nn][AAAA][X]=1/denom;f1ToInf[nn][BBBB][X]=-f2[X]/denom;f1ToInf[nn][CCCC][X]=1/denom;f1ToInf[nn][DDDD][X]=-f1[X]/denom;
	  	f1ToInf[nn][AAAA][Y]=0;      f1ToInf[nn][BBBB][Y]=-f2[Y]/denom;f1ToInf[nn][CCCC][Y]=0;	 f1ToInf[nn][DDDD][Y]=-f1[Y]/denom;
	}
	/**** As right circular cones **************************/
	FOR(i,0,NAP){/**** The template Equator ******************/
		equator[lastLune][0][i][X]=atan(theta[0]) *cos(2*Pi*(float)i/NAP);
		equator[lastLune][0][i][Y]=atan(theta[0]) *sin(2*Pi*(float)i/NAP);
		equator[lastLune][0][i][Z]=1;
	}
	FOR(i,0,NAP){/**** the template Axis  ****************/
		fac=begFac*pow(multiplier,i);
		axes[lastLune][i][Z]=fac;
	}
	FOR(k,0,TH){ /* the template lune *******************/
		FOR(i,0,NA){
			FOR(j,0,NAP){
				fac=begFac*pow(multiplier,i);
				lune[lastLune][k][i][j][X]=fac*atan(theta[k]) *cos(2*Pi*(float)j/NAP);
				lune[lastLune][k][i][j][Y]=fac*atan(theta[k]) *sin(2*Pi*(float)j/NAP);
				lune[lastLune][k][i][j][Z]=fac;
			}
		}
       	}
	FOR(i,0,NR){ /** The template set of rings **/
		FOR(j,0,5){
			FOR(k,0,NAP){
				fac=6*begFac*pow(multiplier,5*i);
				ring[lastLune][0][i][j][k][X]=fac*atan((j+1)*theta[0])*cos(2*Pi*(float)k/NAP);
				ring[lastLune][0][i][j][k][Y]=fac*atan((j+1)*theta[0])*sin(2*Pi*(float)k/NAP);
				ring[lastLune][0][i][j][k][Z]=fac;
			}
		}
	}
	/*******
	For Elliptic & Loxodromic == lune[0].  Klein4 uses lune[1] as well.
	*************/
		FOR(i,0,NA){
			FOR(j,0,NAP){
				pE(lune[0][0][i][j],lune[lastLune][0][i][j],infToF1[0]);
				pE(lune[1][0][i][j],lune[lastLune][0][i][j],infToF1[1]);
			}
		}
		FOR(i,0,NR){
		   FOR(k,0,5){
			FOR(j,0,NAP){
				pE(ring[0][0][i][k][j],ring[lastLune][0][i][k][j],infToF1[0]);
			}
		   }
		}
			FOR(j,0,NAP){
				pE(axes[0][j],axes[lastLune][j],infToF1[0]);
				pE(axes[1][j],axes[lastLune][j],infToF1[1]);
				GETS(frameAxes[0][j],axes[0][j]);
				GETS(frameAxes[1][j],axes[1][j]);
			}
	/***********
	For Klein 4 Group
	*****************/

	/***********
	For Parabolic
	*****************/
	constructTetrahedron();
	newFaces();

}
void drawParaPlane(void){
	int ii, jj;
	FOR(ii,1,30){
			glBegin(GL_LINE_STRIP);
			glVertex3fv(paraPlane[2*ii][0]);
			glVertex3fv(paraPlane[2*ii+1][0]);
			glEnd();
	}
	FOR(jj,1,30){
			glBegin(GL_LINE_STRIP);
			glVertex3fv(paraPlane[0][2*jj]);
			glVertex3fv(paraPlane[0][2*jj+1]);
			glEnd();
	}
}
void drawPlane(void){
	int ii, jj, kk;

	glLineWidth(1);
	glColor3f(.5,.5,.5);
	FOR(ii,0,NLATS){
		glBegin(GL_LINE_STRIP);
		FOR(jj,0,NAP){
			glVertex3fv(fPLats[ii][jj]);
		}
		glEnd();
	}
	FOR(ii,0,NLONGS){
		glBegin(GL_LINE_STRIP);
		FOR(jj,1,NAP){
			if((pow(fPLongs[ii][jj][X],2)+pow(fPLongs[ii][jj][Y],2))>100){
				glEnd();
				while((pow(fPLongs[ii][jj][X],2)+pow(fPLongs[ii][jj][Y],2))>100){
					jj++;
				}
				glBegin(GL_LINE_STRIP);
			}
			glVertex3fv(fPLongs[ii][jj]);
		}
		glEnd();
	}
}
float computeCardinalSpline(float P0,float P1,float P2,float P3,float t){
	float temp[4];

	    /* {-1,  1, -1,  1} {P0} {temp[0]} */
	    /* { 2, -2,  1, -1} {P1} {temp[1]} */
	    /* {-1,  0,  1,  0}*{P2}={temp[2]} */
	    /* { 0,  1,  0,  0} {P3} {temp[3]} */
	    /* Cardinal Matrix K times the x-coordinates of 4 controlPnts  */

	temp[0]= -P0 + P1-P2+P3;
	temp[1]=2*P0-2*P1+P2-P3;
	temp[2]= -P0   +  P2;
	temp[3]=       P1;

	/* {t^3,t^2,t,1} {temp[0]
					  temp[1]
					  temp[2]
					  temp[3]} */

	tempTan=3*t*t*temp[0]+2*t*temp[1]+temp[2];
	tempNorm=6*t*temp[0]+2*temp[1];

	return t*t*t*temp[0]+t*t*temp[1]+t*temp[2]+temp[3];

}
void initBBControlPnts(int chN){
	int ii,jj;

	FOR(jj,0,3){
		controlPnts[0][jj]=aCent[bbIndices[chN][1]][jj];
		controlPnts[lastBB[chN]+1][jj]=aCent[ bbIndices[chN][lastBB[chN]] ][jj];
	}
	FOR(ii,1,lastBB[chN]+1){
		FOR(jj,0,3){controlPnts[ii][jj]=aCent[bbIndices[chN][ii]][jj];}
	}
}
/******************************
void computeFrame(void){
	float magT, magB;
	int ii, jj, kk, mm;
	FOR(mm,1,lastChain+1){
		FOR(ii,0,lastBB[mm]){
			FOR(jj,0,spineKnots){
				CROSS(NormA[mm][ii*spineKnots+jj],TanA[mm][ii*spineKnots+jj],Frame[BIN][mm][ii*spineKnots+jj]);
				magT=NRM(TanA[mm][ii*spineKnots+jj]);
				magB=NRM(Frame[BIN][mm][ii*spineKnots+jj]);
				FOR(kk,0,3){
					Frame[TAN][mm][ii*spineKnots+jj][kk]=TanA[mm][ii*spineKnots+jj][kk]/magT;
					Frame[BIN][mm][ii*spineKnots+jj][kk]/=magB;
				}
				CROSS(Frame[TAN][mm][ii*spineKnots+jj],Frame[BIN][mm][ii*spineKnots+jj],Frame[NOR][mm][ii*spineKnots+jj]);
			}
		}
	}
}
********************************************/
void computeBackBone(void){
	float tParam;
	int ii, jj, kk, mm;
   FOR(mm,1,lastChain+1){
	initBBControlPnts(mm);
	FOR(ii,0,lastBB[mm]){
		FOR(jj,0,spineKnots){
			tParam=((float)jj)/((float)spineKnots);
		FOR(kk,0,3){
	backbone[mm][ii*spineKnots+jj][kk]=computeCardinalSpline(controlPnts[ii][kk],controlPnts[ii+1][kk],controlPnts[ii+2][kk],controlPnts[ii+3][kk],tParam);

	/** insert frame computations here?? **
		TanA[mm][ii*spineKnots+jj][kk]=tempTan;
		NormA[mm][ii*spineKnots+jj][kk]=tempNorm;
	******************************/
			}
		}
	}
   }
}
void constructPrototypeAtomSphere(void){
	int jj, kk;
	float ttheta, pphi;

	FOR(jj,1,NA2+1){
		FOR(kk,0,NA2+1+1){
			ttheta=2*Pi*(float)kk/NA2;
			pphi=Pi*(float)jj/NA2;
			aSphere[jj][kk][X]=cos(Pi/2-pphi)*cos(ttheta);
			aSphere[jj][kk][Y]=cos(Pi/2-pphi)*sin(ttheta);
			aSphere[jj][kk][Z]=sin(Pi/2-pphi);
		}
	}

}
void reConstructMarker(void){
	int ii, jj, kk;
	float rrad=.5, ggap=1.3, sshift;
	float ttheta, pphi;
	ii=0;
}
void constructLox(void){
	int i, j, k, l, d, a, o, r, scN;
	float t1[3],t2[3],fac,begFac,multiplier,temp[3],xr,yi,mult2,fac2;
	multiplier=multiplierK;test1[0]=multiplier;
	begFac=pow(1/multiplier,NAP/2-1);
	fac2=.1;
	mult2=pow(multiplier,fac2);
	FOR(i,0,NA*10){
		fac=begFac*pow(mult2,i);
		t1[X]=fac*atan(theta[0]) *cos(0.f);
		t1[Y]=fac*atan(theta[0]) *sin(0.f);
		t1[Z]=0;
		t2[X]=cos(i*twist*fac2);t2[Y]=sin(i*twist*fac2);t2[Z]=0;
		mul(temp,t1,t2);
		temp[Z]=0.;
		pE(loxCurvePoints[0][i],temp,infToF1[0]);
		GETS(frameLoxCurvePoints[0][i],loxCurvePoints[0][i]);
		temp[Z]=fac;
		pE(loxCurvePoints[1][i],temp,infToF1[0]);
		GETS(frameLoxCurvePoints[1][i],loxCurvePoints[1][i]);
	}
}
void plane(void){
	float rad,cent[2],tempTheta,ttt,LL=5,temp[3]={0,0,0};
	int i, j, k;
	int dif;
	dif=(int)((float)NAP/2);

	/** plane for parabolic case **/
	FOR(i,1,30){
			paraPlane[2*i][0][X]=(-15+i)+cos(Pi/3)*(-15);
			paraPlane[2*i][0][Y]=(-15)*sin(Pi/3);
			paraPlane[2*i+1][0][X]=(-15+i)+cos(Pi/3)*(15);
			paraPlane[2*i+1][0][Y]=(15)*sin(Pi/3);
	}
	FOR(j,1,30){
			paraPlane[0][2*j][X]=(-15+j)-cos(Pi/3)*(-15);
			paraPlane[0][2*j][Y]=(-15)*sin(Pi/3);
			paraPlane[0][2*j+1][X]=(-15+j)-cos(Pi/3)*(15);
			paraPlane[0][2*j+1][Y]=(15)*sin(Pi/3);
	}

	FOR(i,0,NH){	/*** hyp Pencil Circles */
		cent[X]=(fixedPnts[0][0][X]+fixedPnts[0][1][X])/2;
		cent[Y]=(float)i-(float)NH/2;
		rad=pow(pow(cent[X]-fixedPnts[0][0][X],2.f)+
			pow(cent[Y]-fixedPnts[0][0][Y],2.f),.5f);
		if(cent[Y]<=0) {tempTheta=atan(-(fixedPnts[0][0][X]-cent[X])/cent[Y]);}
		else {tempTheta=atan((fixedPnts[0][0][X]-cent[X])/cent[Y]);}
		FOR(j,0,NHAP){
		   if(cent[Y]<=0){/********** in lower half plane **********/
			if(j<=(dif)){/*** first half of arc ****/
				ttt=(Pi/2)+tempTheta*(float)j/(dif);
			}
			else {
				ttt=Pi/2 + tempTheta +(Pi-tempTheta)*(float)(j-dif)/(NHAP-dif);
			}
		   }
		   else {
			if(j<=(dif)){
				ttt=(-Pi/2)-tempTheta*(float)j/(dif);
			}
			else {
				ttt=-Pi/2 - tempTheta -(Pi-tempTheta)*(float)(j-dif)/(NHAP-dif);
			}
		   }	/** 1/2- way model of pencil of hyperbolic circles in gen position **/
		   if(cent[Y]<=0){
			if(j<=(dif)){
				ttt=(Pi/2-tempTheta)+2*tempTheta*(float)j/(dif);
			}
			else {
				ttt=Pi/2 + tempTheta +(2*Pi-2*tempTheta)*(float)(j-dif)/(NHAP-dif);
			}
		   }
		   else {
			if(j<=(dif)){
				ttt=(-Pi/2+tempTheta)-2*tempTheta*(float)j/(dif);
			}
			else {
				ttt=-Pi/2-tempTheta-(2*Pi-2*tempTheta)*(float)(j-dif)/(NHAP-dif);
			}
		   }	/** pencil of hyperbolic circles in gen position **/
			hypPencil[i][j][X]=rad*cos(ttt) + cent[X];
			hypPencil[i][j][Y]=rad*sin(ttt) + cent[Y];
			hypPencil[i][j][Z]=0;
		}
	}

	FOR(i,0,NLATS){
	       FOR(j,0,NAP){
			sLats[i][j][X]=cos((float)j/(NAP-1)*2*Pi)*sin((float)(i+1)/(NLATS+1)*Pi);
			sLats[i][j][Y]=sin((float)j/(NAP-1)*2*Pi)*sin((float)(i+1)/(NLATS+1)*Pi);
			sLats[i][j][Z]=cos((float)(i+1)/(NLATS+1)*Pi);
		}
	}
	FOR(i,0,NLONGS){
	        FOR(j,0,NAP){
			sLongs[i][j][X]=cos((float)i/(NLONGS)*2*Pi)*sin((float)(j)/(NAP-1)*Pi);
			sLongs[i][j][Y]=sin((float)i/(NLONGS)*2*Pi)*sin((float)(j)/(NAP-1)*Pi);
			sLongs[i][j][Z]=cos((float)(j)/(NAP-1)*Pi);
		}
	}
	FOR(i,0,NLATS){
	        FOR(j,0,NAP){
			inverseStereoProjection(temp,sLats[i][j]);
			pE(pLats[i][j],temp,infToF1[0]);
		}
	}
	FOR(i,0,NLONGS){
	        FOR(j,0,NAP){
			inverseStereoProjection(temp,sLongs[i][j]);
			/* hack to eliminate point mapping to infinity */
			if(temp[X]!=1 || temp[Y]!=0) pE(pLongs[i][j],temp,infToF1[0]);
			else GETS(pLongs[i][j],pLongs[i][j-1]);
		}
	}
	FOR(i,0,NLATS){
	        FOR(j,0,NAP){
			GETS(fPLats[i][j],pLats[i][j]);
		}
	}
	FOR(i,0,NLONGS){
	        FOR(j,0,NAP){
			GETS(fPLongs[i][j],pLongs[i][j]);
		}
	}

}
void initializeFigures(){
	int ii, jj, kk, aaa, ppp, rr, c, t,f;
		FOR(ii,0,NA*10){
			GETS(frameLoxCurvePoints[0][ii],loxCurvePoints[0][ii]);
			GETS(frameLoxCurvePoints[1][ii],loxCurvePoints[1][ii]);
		}
		FOR(aaa,0,NA){/**** single layer only ***/
			FOR(ppp,0,NAP){
				GETS(frameLune[0][0][aaa][ppp],lune[0][0][aaa][ppp]);
			}
		}

/**************************************/
	FOR(ii,0,NLATS){
	    FOR(jj,0,NAP){
			GETS(fPLats[ii][jj],pLats[ii][jj]);
		}
	}
	FOR(ii,0,NLONGS){
	    FOR(jj,0,NAP){
			GETS(fPLongs[ii][jj],pLongs[ii][jj]);
		}
	}
		FOR(ppp,0,NAP){
			GETS(frameAxes[0][ppp],axes[0][ppp]);
		}
	FOR(rr,0,NR){
		FOR(c,0,5){
	       	FOR(ppp,0,NAP){
				GETS(frameRing[0][rr][c][ppp],ring[0][0][rr][c][ppp]);
			}
		}
	}
   	FOR(t,0,2){
	   FOR(f,0,4){
		FOR(ii,0,3){
			FOR(jj,0,3){
				FOR(kk,0,NAP){
					GETS(frameTetra[t][f][ii][jj][kk],tetra[t][f][ii][jj][kk]);
				}
			}
		}
	   }
	}

}
/**********************************************************************/
/**********************************************************************/
/****** Miscelanous functions - do various things, not calc or draw****/
/**********************************************************************/
/**********************************************************************/
/**********************************************************************/
/************** Spherical Geometry subroutines ************************/
/**********************************************************************/
/****************/
void crossProduct(float res[3],float p1[3],float p2[3]){
	res[X]=p1[Y]*p2[Z]-p1[Z]*p2[Y];
	res[Y]=-p1[X]*p2[Z]+p1[Z]*p2[X];
	res[Z]=p1[X]*p2[Y]-p1[Y]*p2[X];
}
void rotateAboutAxis(float res[3],float axis[3],float pnt[3],float theta10){
	int ii, jj, kk;
	float aXp[3], AdP, nrm, tmp[3], ax[3], projCoeff;

	/**normalize axis**/nrm=NRM(axis);FOR(ii,0,3){ax[ii]=axis[ii]/nrm;}
	projCoeff=DOT(ax,pnt);/*==DOT(axis,pnt)/(nrm*nrm);*/
	crossProduct(aXp,ax,pnt);
	/**normalize aXp**/  nrm=NRM(aXp); FOR(ii,0,3){aXp[ii]/=nrm;}
	/**coeff for aXp. Needs to be same length as P==(pnt[ii]-AdP*ax[ii]),
	so that the aXp term lies on the same circle perp to and centered on axis
	as does P==(pnt[ii]-AdP*ax[ii])**/
	FOR(ii,0,3){tmp[ii]=pnt[ii]-projCoeff*ax[ii];}nrm=NRM(tmp);
	FOR(ii,0,3){
  res[ii]=projCoeff*ax[ii]+sin(theta10)*nrm*aXp[ii]+cos(theta10)*(pnt[ii]-projCoeff*ax[ii]);
	}
}
void multMatrices3D(float res[3][3],float preMat[3][3],float postMat[3][3]){
	/******************************************
	res[0][0]=preMat[0][0]*postMat[0][0]+preMat[0][1]*postMat[1][0]+preMat[0][2]*postMat[2][0];
	res[0][1]=preMat[0][0]*postMat[0][1]+preMat[0][1]*postMat[1][1]+preMat[0][2]*postMat[2][1];
	res[0][2]=preMat[0][0]*postMat[0][2]+preMat[0][1]*postMat[1][2]+preMat[0][2]*postMat[2][2];

	res[1][0]=preMat[1][0]*postMat[0][0]+preMat[1][1]*postMat[1][0]+preMat[1][2]*postMat[2][0];
	res[1][1]=preMat[1][0]*postMat[0][1]+preMat[1][1]*postMat[1][1]+preMat[1][2]*postMat[2][1];
	res[1][2]=preMat[1][0]*postMat[0][2]+preMat[1][1]*postMat[1][2]+preMat[1][2]*postMat[2][2];

	res[2][0]=preMat[2][0]*postMat[0][0]+preMat[2][1]*postMat[1][0]+preMat[2][2]*postMat[2][0];
	res[2][1]=preMat[2][0]*postMat[0][1]+preMat[2][1]*postMat[1][1]+preMat[2][2]*postMat[2][1];
	res[2][2]=preMat[2][0]*postMat[0][2]+preMat[2][1]*postMat[1][2]+preMat[2][2]*postMat[2][2];
	******************************************/

	FOR(ii,0,3){
		FOR(jj,0,3){
			res[ii][jj]=preMat[ii][0]*postMat[0][jj]+preMat[ii][1]*postMat[1][jj]+preMat[ii][2]*postMat[2][jj];
		}
	}
}
void multMatrices4D(float res[4][4],float preMat[4][4],float postMat[4][4]){
	FOR(ii,0,4){
		FOR(jj,0,4){
			res[ii][jj]=preMat[ii][0]*postMat[0][jj]+preMat[ii][1]*postMat[1][jj]+preMat[ii][2]*postMat[2][jj]+preMat[ii][3]*postMat[3][jj];
		}
	}
}
/***********************************************************************/
/**********************************************************************/
/* invert - Inverts a aff matrix                                */
/**********************************************************************/
void invert(float *inv, float *mat){
 inv[0] = mat[0];
 inv[1] = mat[4];
 inv[2] = mat[8];
 inv[3] = mat[3];
 inv[4] = mat[1];
 inv[5] = mat[5];
 inv[6] = mat[9];
 inv[7] = mat[7];
 inv[8] = mat[2];
 inv[9] = mat[6];
 inv[10] = mat[10];
 inv[11] = mat[11];
 inv[12] = -mat[0]*mat[12] - mat[1]*mat[13] - mat[2]*mat[14];
 inv[13] = -mat[4]*mat[12] - mat[5]*mat[13] - mat[6]*mat[14];
 inv[14] = -mat[8]*mat[12] - mat[9]*mat[13] - mat[10]*mat[14];
 inv[15] = mat[15];

}
/**********************************************************************/
/* updatecontrail - This gets the places to make the contrail.        */
/**********************************************************************/
void updatecontrail(void){
  static int flightpt=0;
  float temp[16];
  if(flightpt>=LONGFLIGHT){
      flightpt=0;
    }
    glGetFloatv(GL_MODELVIEW_MATRIX,temp);
    invert(contrail[flightpt], temp);
    flightpt++;
}
/**********************************************************************/
/* odometer - computes the framerate                                  */
/**********************************************************************/
#ifdef IRIX
float odometer(void){ /* right now we're not using the odometer */
    static struct timeval tv0;  /* time when program started */
    struct timeval tv;
    gettimeofday(&tv, NULL);
    if(tv0.tv_sec == 0) tv0 = tv;  /* initialize tv0 */
    return (tv.tv_sec - tv0.tv_sec) + (tv.tv_usec - tv0.tv_usec)*1e-6;
}
#endif
#ifdef OSXMac
float odometer(void){ /* right now we're not using the odometer */
    static struct timeval tv0;  /* time when program started */
    struct timeval tv;
    gettimeofday(&tv, NULL);
    if(tv0.tv_sec == 0) tv0 = tv;  /* initialize tv0 */
    return (tv.tv_sec - tv0.tv_sec) + (tv.tv_usec - tv0.tv_usec)*1e-6;
}
#endif
/**********************************************************************/
/**********************************************************************/

/** DNA stuff **/
#define bgntmesh() 	glBegin(GL_TRIANGLE_STRIP)
#define endtmesh()	glEnd()
#define v3f		glVertex3fv
#define bgnline()	glBegin(GL_LINES)
#define bgnlineS()	glBegin(GL_LINE_STRIP)
#define endline()	glEnd()
#define bgnpoint()	glBegin(GL_POINTS)
#define endpoint()	glEnd()

/**********************************************************************/
/**********************************************************************/
/******** Drawing functions *******************************************/
/**********************************************************************/
/**********************************************************************/
/**********************************************************************/
/**********************************************************************/
/**********************************************************************/
/* drawcube() - Draws a cube for visual reference                     */
/**********************************************************************/
void drawcube(void){
      glColor3f(0.0,0.0,1.0);
      glLineWidth(3);
      glutWireCube(2);
      /*glutWireSphere(1,10,10); */
}
void drawsphere(void){
      glColor3f(0.0,0.0,1.0);
      /**glLineWidth(3);
      glutWireSphere(.1,10,10); **/

          glTranslatef(1,10 ,.25);
          glScalef(1, (float)20, 1);
      glutWireSphere(.1,10,10);

}
/**********************************************************************/
/* drawcontrail - Draws the contrail                                  */
/**********************************************************************/
void drawcontrail(void){
   int ii;
   glLineWidth(1);
   for(ii=0; ii<=LONGFLIGHT; ii++){
     glPushMatrix();
     glMultMatrixf(contrail[ii]);
       glBegin(GL_LINES);

         glColor3f(1,0,0);
         glVertex3f( 1, 0, 0);
         glVertex3f(-1, 0, 0);

         glColor3f(0,1,0);
         glVertex3f( 0, 1, 0);
         glVertex3f( 0,-1, 0);

         glColor3f(0,0,1);
         glVertex3f( 0, 0, 1);
         glVertex3f( 0, 0,-1);
       glEnd();
    glPopMatrix();
    }
}
/**********************************************************************/
void Hlight2(int inx){
	float spec; int colorT;
	float amb= .33, pwr= 32.; int  brt=255;

	switch(backBone){
		case LASTBBONE:
		case 0:
		case 4:/*backbone core only*/
		case 1: colorT=atomInfoNum[inx][ATOTYP]; /**/break;/**/
		/* color by hydrophobicity */
        case 2: /** connectivity only **/
 		case 3: /** connectivity and spheres */
       	if(atomInfoNum[inx][ATOTYP]!=0){
			   switch(atomInfoNum[inx][RESTYN]){ /* set color for hydrophilic/phobic side chain */
				case ALA: /* hydrophobic */
				case CYS:
				case ILE:
				case LEU:
				case MET:
				case PHE:
				case PRO:
				case TRP:
				case TYR:
				case VAL:
		   			colorT=PHOB; break;
				case ARG:
				case HIS:
				case LYS:
       				colorT=CHRGP; break;
				case ASP:  /* charged */
				case GLU:
				 	colorT=CHRGN; break;
				case ASN:
				case GLN:
				case GLY:
				case SER:
				case THR:
       				colorT=PHIL; break;
	   		   }/*end switch(residueTyp)*/
			}
			else{ colorT=0;}/*this for marker sphere*/	break;
	}/* end switch backBone */
		rred=col[colorT][0];
		ggre=col[colorT][1];
		bblu=col[colorT][2];
    /*note that lmb is computed in drawSphere()*/
	/*****************/
	lmb = MAX(lmb,amb);
    spec=brt*( 1. - pwr + pwr*lmb);
    rred = MAX(lmb*rred, spec);
    ggre = MAX(lmb*ggre, spec);
    bblu = MAX(lmb*bblu, spec);
    /***********************/

    if(opaque) opaq=.7;  else opaq=1.0;

    glColor4f(rred,ggre,bblu,opaq);
}
void drawSphere(int nn){
	int t = atomInfoNum[nn][ATOTYP], ii, jj;
	float r = atomRadius[ t ];
	float atompos[3];

	GETS( atompos, aCent[ nn ] );
	for(ii = 0; ii < NA2; ii += 2) {
	    glBegin( GL_TRIANGLE_STRIP );
	    for(jj = 1; jj < NA2-1; jj++) {
			float *v1 = &aSphere[jj][ii][0];  /* so v1[0]=x, v1[1]=y, etc.*/
			float *a1 = &aSphere[jj][ii+1][0];

			lmb = fabs( luxx[0]*v1[0] + luxx[1]*v1[1] + luxx[2]*v1[2] );

			Hlight2( nn );

			glVertex3f( atompos[0] + r * v1[0],
			    atompos[1] + r * v1[1],
			    atompos[2] + r * v1[2] );

			glVertex3f( atompos[0] + r * a1[0],
			    atompos[1] + r * a1[1],
			    atompos[2] + r * a1[2] );
	    }
	    glEnd();
	}
}
void drawCoreOnly(void){
	int ii, jj;
	FOR(ii,1,lastChain+1){
		FOR(jj,1,lastBB[ii]+1){
			drawSphere(bbIndices[ii][jj]);
		}
	}
}
void drawMarkerSphere(void){
	int ii, jj;
	float r=atomRadius[0];
	float atompos[3];
	if(bigMarker){r=.5*10;}

	GETS( atompos, aCent[ alphaCIndices[targetAC] ] );
	for(ii = 1; ii < NA2+1; ii += 2) {
	    glBegin( GL_TRIANGLE_STRIP );
	    for(jj = 1; jj < NA2-1; jj++) {
			float *v1 = &aSphere[jj][ii][0];  /* so v1[0]=x, v1[1]=y, etc.*/
			float *a1 = &aSphere[jj][ii+1][0];

			lmb = fabs( luxx[0]*v1[0] + luxx[1]*v1[1] + luxx[2]*v1[2] );

			Hlight2( 0 );

			glVertex3f( atompos[0] + r * v1[0],
			    atompos[1] + r * v1[1],
			    atompos[2] + r * v1[2] );

			glVertex3f( atompos[0] + r * a1[0],
			    atompos[1] + r * a1[1],
			    atompos[2] + r * a1[2] );
	    }
	    glEnd();
	}
}
/**********************************************************************/
/* drawall - draws the main viewport.  this is where one should toggle*/
/*    things we want to draw.                                         */
/**********************************************************************/
void drawTetras(){
	int t, f, i, j, k;
   FOR(t,0,2){/*frameTetra[t][f][i][j][k]*/
     	FOR(f,0,4){
		FOR(i,0,3){
			FOR(j,0,2){
				switch(i){
					case 0:
					if(j==0) {
						switch(f){
							case 0:glColor4f(0.0,1.0,1.0,0.7);break;
							case 1:glColor4f(0.0,0.5,1.0,0.7);break;
							case 2:glColor4f(1.0,1.0,0.0,0.7);break;
							case 3:glColor4f(1.0,0.5,0.0,0.7);break;
						}/* f,0 */
					}
					else {
						switch(f){
							case 0:glColor4f(1.0,1.0,0.0,0.7);break;
							case 1:glColor4f(1.0,0.5,0.0,0.7);break;
							case 2:glColor4f(0.0,1.0,1.0,0.7);break;
							case 3:glColor4f(1.0,0.7,0.7,0.7);break;
						}/* f,1 */
					}
					break;

					case 1:
					if(j==1){
						switch(f){
							case 0:glColor4f(0.0,1.0,1.0,0.7);break;
							case 1:glColor4f(0.0,0.5,1.0,0.7);break;
							case 2:glColor4f(1.0,1.0,0.0,0.7);break;
							case 3:glColor4f(1.0,0.5,0.0,0.7);break;
						}/* f,0 */
					}
					else {
						switch(f){
							case 0:glColor4f(1.0,0.3,0.3,0.7);break;
							case 1:glColor4f(1.0,0.5,0.5,0.7);break;
							case 2:glColor4f(1.0,1.3,0.3,0.7);break;
							case 3:glColor4f(0.0,0.5,1.0,0.7);break;
						}/* f,2 */
					}
					break;

					case 2:
					if(j==1){
						switch(f){
							case 0:glColor4f(1.0,0.3,0.3,0.7);break;
							case 1:glColor4f(1.0,0.5,0.5,0.7);break;
							case 2:glColor4f(1.0,1.3,0.3,0.7);break;
							case 3:glColor4f(0.0,0.5,1.0,0.7);break;
						}/* f,2 */
					}
					else {
						switch(f){
							case 0:glColor4f(1.0,1.0,0.0,0.7);break;
							case 1:glColor4f(1.0,0.5,0.0,0.7);break;
							case 2:glColor4f(0.0,1.0,1.0,0.7);break;
							case 3:glColor4f(1.0,0.7,0.7,0.7);break;
						}/* f,1 */
					}
					break;
				}
				if((band==1 && f!=1)) glBegin(GL_TRIANGLE_STRIP);
				FOR(k,0,NAP){
					/*for dragon teeth*/ /**/if((k-(int)((float)k/4)*4==0 && band==1) && (f!=1))/**/
					glVertex3fv(frameTetra[t][f][i][j][k]);
					glVertex3fv(frameTetra[t][f][i][j+1][k]);
					if((k-(int)((float)k/4)*4==2 || k==NAP-1) /**/&& band==0/**/) glEnd();
				}
				if(band==1) glEnd();
			}
		}
     	}/**end FOR(f,0,numFaces)...**/
   	FOR(f,0,4){
		FOR(i,0,3){
			FOR(j,0,3){
				glBegin(GL_LINES);
				bgnline();
				FOR(k,0,NAP){
					glVertex3fv(frameTetra[t][f][i][j][k]);
				}
				glEnd();
			}
		}
	}
  }
}




/*******************
FOR(mm,1,lastChain+1){
	initBBControlPnts(mm);
	FOR(ii,0,lastBB[mm]){
		FOR(jj,0,spineKnots){
			tParam=((float)jj)/((float)spineKnots);
		FOR(kk,0,3){
	backbone[mm][ii*spineKnots+jj][kk]=computeCardinalSpline(controlPnts[ii][kk],controlPnts[ii+1][kk],controlPnts[ii+2][kk],controlPnts[ii+3][kk],tParam);

		TanA[mm][ii*spineKnots+jj][kk]=tempTan;
		NormA[mm][ii*spineKnots+jj][kk]=tempNorm;

			}
		}
	}
   }
***************************/
   /****************************
   Need to declare control points for each chain!!!
   and run initBBControlPnts(mm); just once for all chains....
   No, have already declared indices for backbone atoms
   w/ chain# and position# on chain specified.
   **********************************
computeFrameSinglePnt(int ch, int cp, int sk){
	int kk;
	float tempT[3], tempN[3], tempB[3], tempBB[3], tParam, magT, magN;
	tParam=(float)sk/((float)spineKnots)
	FOR(kk,0,3){
		computeCardinalSpline2(controlPnts[ch][cp][sk][kk],controlPnts[ch][cp+1][sk][kk],controlPnts[ch][cp+2][sk][kk],controlPnts[ch][cp+3][sk][kk],tParam);
		tempT[kk]=tempTan;
		tempN[kk]=tempNorm;
		tempBB[kk]=tempBB;
	}
	magT=NRM(tempT);
	magN=NRM(tempN);
	FOR(kk,0,3){
		tempT[kk]/=magT;
		tempN[kk]/=magN;
	}
	CROSS(tempT,tempN,tempB);
	CROSS(tempB,tempT,tempN);
}
************************/
void drawSingleFrame(int chN,int atN,int jtN){
	float tempT[3], tempN[3], tempB[3], tempPnt[3], tempTip[3];
	float magT, magN, tParam, p0, p1, p2, p3;
	int xyz;/*ch(ain)N(umber), at(om)N(on each chain), j(oin)tN(btwn atoms)*/
				tParam=(float)jtN/((float)spineKnots);
				FOR(xyz,0,3){
					if(atN==0){
						p0=p1=aCent[bbIndices[chN][atN]][xyz];
						p2=aCent[bbIndices[chN][atN+1]][xyz];
						p3=aCent[bbIndices[chN][atN+2]][xyz];
					}
					if(atN==lastBB[chN]-2){
						p0=aCent[bbIndices[chN][atN-1]][xyz];
						p1=aCent[bbIndices[chN][atN]][xyz];
						p2=p3=aCent[bbIndices[chN][atN+1]][xyz];
					}
					if(atN>0 && atN<lastBB[chN]-2){
						p0=aCent[bbIndices[chN][atN-1]][xyz];
						p1=aCent[bbIndices[chN][atN]][xyz];
						p2=aCent[bbIndices[chN][atN+1]][xyz];
						p3=aCent[bbIndices[chN][atN+2]][xyz];
					}
					computeCardinalSpline2(p0,p1,p2,p3,tParam);
					tempT[xyz]=tempTan;
					tempN[xyz]=tempNorm;
					tempPnt[xyz]=tempBB;
				}
				magT=NRM(tempT);
				magN=NRM(tempN);
				FOR(xyz,0,3){
					tempT[xyz]/=magT;
					tempN[xyz]/=magN;
				}
				CROSS(tempT,tempN,tempB);
				CROSS(tempB,tempT,tempN);

				/** compute point at tip of tan component **/
				FOR(xyz,0,3){
					tempTip[xyz]=tempPnt[xyz]+tempT[xyz];
				}
				glColor3f(1.0,0.0,0.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(tempPnt);
		    	glVertex3fv(tempTip);
	    		glEnd();
				/** compute point at tip of norm component **/
				FOR(xyz,0,3){
					tempTip[xyz]=tempPnt[xyz]+tempN[xyz];
				}
				glColor3f(0.0,1.0,0.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(tempPnt);
		    	glVertex3fv(tempTip);
	    		glEnd();
				/** compute point at tip of binorm component **/
				FOR(xyz,0,3){
					tempTip[xyz]=tempPnt[xyz]+tempB[xyz];
				}
				glColor3f(1.0,0.0,1.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(tempPnt);
		    	glVertex3fv(tempTip);
	    		glEnd();
}
void computeCardinalSpline2(float P0,float P1,float P2,float P3,float t){
	float temp[4];

	    /* {-1,  1, -1,  1} {P0} {temp[0]} */
	    /* { 2, -2,  1, -1} {P1} {temp[1]} */
	    /* {-1,  0,  1,  0}*{P2}={temp[2]} */
	    /* { 0,  1,  0,  0} {P3} {temp[3]} */
	    /* Cardinal Matrix K times the x-coordinates of 4 controlPnts  */

	temp[0]= -P0 + P1-P2+P3;
	temp[1]=2*P0-2*P1+P2-P3;
	temp[2]= -P0   +  P2;
	temp[3]=       P1;

	/* {t^3,t^2,t,1} {temp[0]
					  temp[1]
					  temp[2]
					  temp[3]} */

	tempTan=3*t*t*temp[0]+2*t*temp[1]+temp[2];
	tempNorm=6*t*temp[0]+2*temp[1];
	tempBB=t*t*t*temp[0]+t*t*temp[1]+t*temp[2]+temp[3];

}
void drawFrame(void){
	float tempT[3], tempN[3], tempB[3], tempPnt[3], tempTip[3];
	float magT, magN, tParam, p0, p1, p2, p3;
	int atN, jtN, xyz, chN;/*ch(ain)N(umber), at(om)N(on each chain), j(oin)tN(btwn atoms)*/
	FOR(chN,1,lastChain+1){
		FOR(atN,0,lastBB[chN]-1){
			FOR(jtN,0,spineKnots){
				tParam=(float)jtN/((float)spineKnots);
				FOR(xyz,0,3){
					if(atN==0){
						p0=p1=aCent[bbIndices[chN][atN]][xyz];
						p2=aCent[bbIndices[chN][atN+1]][xyz];
						p3=aCent[bbIndices[chN][atN+2]][xyz];
					}
					if(atN==lastBB[chN]-2){
						p0=aCent[bbIndices[chN][atN-1]][xyz];
						p1=aCent[bbIndices[chN][atN]][xyz];
						p2=p3=aCent[bbIndices[chN][atN+1]][xyz];
					}
					if(atN>0 && atN<lastBB[chN]-2){
						p0=aCent[bbIndices[chN][atN-1]][xyz];
						p1=aCent[bbIndices[chN][atN]][xyz];
						p2=aCent[bbIndices[chN][atN+1]][xyz];
						p3=aCent[bbIndices[chN][atN+2]][xyz];
					}
					computeCardinalSpline2(p0,p1,p2,p3,tParam);
					tempT[xyz]=tempTan;
					tempN[xyz]=tempNorm;
					tempPnt[xyz]=tempBB;
				}
				magT=NRM(tempT);
				magN=NRM(tempN);
				FOR(xyz,0,3){
					tempT[xyz]/=magT;
					tempN[xyz]/=magN;
				}
				CROSS(tempT,tempN,tempB);
				CROSS(tempB,tempT,tempN);

				/** compute point at tip of tan component **/
				FOR(xyz,0,3){
					tempTip[xyz]=tempPnt[xyz]+tempT[xyz];
				}
				glColor3f(1.0,0.0,0.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(tempPnt);
		    	glVertex3fv(tempTip);
	    		glEnd();
				/** compute point at tip of norm component **/
				FOR(xyz,0,3){
					tempTip[xyz]=tempPnt[xyz]+tempN[xyz];
				}
				glColor3f(0.0,1.0,0.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(tempPnt);
		    	glVertex3fv(tempTip);
	    		glEnd();
				/** compute point at tip of binorm component **/
				FOR(xyz,0,3){
					tempTip[xyz]=tempPnt[xyz]+tempB[xyz];
				}
				glColor3f(1.0,0.0,1.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(tempPnt);
		    	glVertex3fv(tempTip);
	    		glEnd();
			}
		}
	}
}
/*************************
void drawFrame2(void){
	float tempnt[3];
	int ii, jj, kk, mm;
	FOR(mm,1,lastChain+1){
		FOR(ii,0,lastBB[mm]){
			jj=0;/**FOR(jj,0,spineKnots){**/
				/** draw each frame component (if only for jj=0, then only at beginning of segment)**/
				/** compute point at tip of tan component **
				FOR(kk,0,3){
					tempnt[kk]=backbone[mm][ii*spineKnots+jj][kk]+Frame[TAN][mm][ii*spineKnots+jj][kk];
				}
				glColor3f(1.0,0.0,0.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(backbone[mm][ii*spineKnots+jj]);
		    	glVertex3fv(tempnt);
	    		glEnd();
				/** compute point at tip of norm component **
				FOR(kk,0,3){
					tempnt[kk]=backbone[mm][ii*spineKnots+jj][kk]+Frame[NOR][mm][ii*spineKnots+jj][kk];
				}
				glColor3f(0.0,1.0,0.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(backbone[mm][ii*spineKnots+jj]);
		    	glVertex3fv(tempnt);
	    		glEnd();
				/** compute point at tip of norm component **
				FOR(kk,0,3){
					tempnt[kk]=backbone[mm][ii*spineKnots+jj][kk]+Frame[BIN][mm][ii*spineKnots+jj][kk];
				}
				glColor3f(1.0,0.0,1.0);
				glBegin(GL_LINE_STRIP);
		    	glVertex3fv(backbone[mm][ii*spineKnots+jj]);
		    	glVertex3fv(tempnt);
	    		glEnd();
			/**}**
		}
	}
}*******************************/
void drawBackbone(void){
	int ii, mm,nn;
	glColor3f(col[BKBN][0],col[BKBN][1],col[BKBN][2]);
  if(!selectChains){
	  FOR(mm,1,lastChain+1){
        if(backBone==LASTBBONE){glColor3f(col[100+mm][0],col[100+mm][1],col[100+mm][2]);}
	    glBegin(GL_LINE_STRIP);
	    FOR(ii,0,(lastBB[mm]-1)*spineKnots){
		    glVertex3fv(backbone[mm][ii]);
	    }
	    glEnd();
      }
  }
  else{
	  mm=whichChain;
	  if(mm>1){nn=lastBB[mm-1]-1;}
	  else {nn=0;}
	  glColor3f(col[100+mm][0],col[100+mm][1],col[100+mm][2]);
	  	    glBegin(GL_LINE_STRIP);
	  	    FOR(ii,nn,(lastBB[mm]-1)*spineKnots){
	  		    glVertex3fv(backbone[mm][ii]);
	  	    }
	    glEnd();
  }
}
void alphaBeta(int begNum){
	glBegin(GL_LINES);
	glVertex3fv(aCent[begNum+1]);
	glVertex3fv(aCent[begNum+4]);
	glEnd();
}
void carbonOxy(int begNum){
	glBegin(GL_LINES);
	glVertex3fv(aCent[begNum+2]);
	glVertex3fv(aCent[begNum+3]);
	glEnd();
}
void slc(int begNum,int countNum){
	int ii;
	glBegin(GL_LINE_STRIP);
	FOR(ii,begNum,begNum+countNum+1){
		glVertex3fv(aCent[ii]);
	}
	glEnd();
}
void dlc(int begNum,int countNum){
	int ii;
	glBegin(GL_LINES);
	FOR(ii,begNum,begNum+countNum){
		glVertex3fv(aCent[ii]);
		glVertex3fv(aCent[ii+2]);
   }
	glEnd();
}
void drawconnections(int alphaCNum, int residueTyp){
	int bN1, bN2=0, cN1=0, cN2=0, b, colorT;
	b=alphaCNum-1;
	bN1=b+4; bN2=0; cN1=0; cN2=0;
 // cerr << "drawconnections.\n";

			   switch(residueTyp){ /* set color for hydrophilic/phobic side chain */
				case ALA: /* hydrophobic */
				case CYS:
				case ILE:
				case LEU:
				case MET:
				case PHE:
				case PRO:
				case TRP:
				case TYR:
				case VAL:
		   			colorT=PHOB; break;
				case ARG:
				case HIS:
				case LYS:
       				colorT=CHRGP; break;
				case ASP:  /* charged */
				case GLU:
				 	colorT=CHRGN; break;
				case ASN:
				case GLN:
				case GLY:
				case SER:
				case THR:
       				colorT=PHIL; break;
	   		   }/*end switch(residueTyp)*/

 //  cerr << "drawconnections1.\n";
	glColor3f(col[colorT][0],col[colorT][1],col[colorT][2]);

	carbonOxy(b);
 //   cerr << "drawconnections2.\n";
	if(residueTyp!=GLY){/* beginning for all by GLY */
 //  cerr << "drawconnections3.\n";
		alphaBeta(b);
	}
 //  cerr << "drawconnections4.\n";
	switch(residueTyp){
		case GLY:
		case ALA: break;
		case ARG: cN1=5; bN2=b+8; cN2=1; break;
		case ASN:
		case ASP:
		case LEU: cN1=2; bN2=b+5; cN2=1; break;
		case CYS:
		case SER: cN1=1; break;
		case GLN:
		case GLU: cN1=3; bN2=b+6; cN2=1; break;
		case HIS: cN1=2; bN2=b+5; cN2=3; slc(b+8,1); break;
		case ILE: cN1=1; bN2=b+4; cN2=2; break;
		case LYS: cN1=4; break;
		case MET: cN1=3; break;
		case PHE: cN1=2; bN2=b+5; cN2=4; slc(b+9,1); break;
		case PRO: cN1=2;
			glBegin(GL_LINES);
			glVertex3fv(aCent[b]);
			glVertex3fv(aCent[b+6]);
			glEnd();
			break;
		case THR:
		case VAL: cN1=1; bN2=b+4; cN2=1; break;
		case TRP: cN1=2; bN2=b+5; cN2=3; dlc(b+9,3); slc(b+12,1);
			glBegin(GL_LINES);
			glVertex3fv(aCent[b+7]);
			glVertex3fv(aCent[b+10]);
			glVertex3fv(aCent[b+8]);glVertex3fv(aCent[b+9]);
			glEnd();/**/
			break;
		case TYR: cN1=2; bN2=b+5; cN2=4; slc(b+9,2); break;
	}
//   cerr << "drawconnections5.\n";
	/**/slc(bN1,cN1);/**/
	/**/if(cN2!=0){dlc(bN2,cN2);}/**/
 //  cerr << "drawconnectionsEnd.\n";
}
void drawallllll(void){
	switch(chapmat){
		case 0:
			glPushMatrix();
			glMultMatrixf(aff);
			drawall();
			glPopMatrix();
			break;
		case 1:
			glPushMatrix();
			glMultMatrixf(aff2);
			drawall2();
			glPopMatrix();
			break;
		case 2:
		case 3:
			glPushMatrix();
			glMultMatrixf(aff);
			drawall();
			glPopMatrix();

			glPushMatrix();
			glMultMatrixf(aff2);
			drawall2();
			glPopMatrix();
			break;
	}
}
void drawall(void){/*AKG*/
	int ii, jj, kk, ww/*hichLune*/, begLune=0, endLune=1, bandI;

if(axx) drawAxes(); if(changeCent) drawSphere(lastCent);/*mark the in-front-of-you spot  */
	/*plane color*/if(gnd) glColor3f(0,0,0); else {glColor3f(.4,.4,.4);}
	switch(pla){
		case 0:
			break;

		case 1:
			FOR(ii,0,(int)NH/2){ /*** Hyperbolic Pencil of Circles ******/
				glBegin(GL_LINE_STRIP);
				FOR(jj,0,NHAP){
					glVertex3fv(frameHypPencil[ii][jj]);
				}
				glEnd();
			}
			break;

		case 2:
			if(sceneType!=PARABOLIC){drawPlane();}
			else {drawParaPlane();}
			break;
	}
	if(sceneType==LOXODROMIC){ /* loxodromes on lune and in plane */
			glLineWidth(3); /*thick*/
			loxPlaneColor;
			glBegin(GL_LINES);
			FOR(ii,0,loxLength){
				glVertex3fv(frameLoxCurvePoints[0][ii]);
			}
			glEnd();
			glBegin(GL_LINES);
			FOR(ii,0,loxLength){
				glVertex3fv(frameLoxCurvePoints[1][ii]);
			}
			glEnd();
	}
	if(sceneType==PARABOLIC){
			drawTetras();
	}
	if(sceneType==MOLECULE){
  //   cerr << "molecule1.\n";
				/**if(!mov) drawSingleFrame(movieCounter);**/
				if(fram1) drawFrame();
   //  cerr << "molecule2.\n";
			if(backBone!=showHydrophob){FOR(ii,1,lastAtom+1){drawSphere(ii);};}
			else{if(backBone==showHydrophob){drawCoreOnly();}}
			glLineWidth(3);
   //  cerr << "molecule3.\n";
			switch(backBone){
				case 0: /*show nothing but atoms*/ break;
				case 1: 
   //  cerr << "molecule3.1.\n";
          drawBackbone();/** draw peptide backbone and marker connections**/
   // cerr << "molecule3.1.1.\n";
					drawMarkerSphere();
   // cerr << "molecule3.1.2.\n";
					if(!noConnections){
   // cerr << "molecule3.1.3.\n";
						drawconnections(alphaCIndices[targetAC],atomInfoNum[ alphaCIndices[targetAC] ][RESTYN]);/**/
					}
					break;
				case 2:/*display hydrophobicity for all atoms, w/o spheres then with*/
   // cerr << "molecule3.2.\n";
				case 3:
   // cerr << "molecule3.3.\n";
          drawBackbone();
					drawMarkerSphere();
					if(!noConnections){
       					FOR(ii,1,lastResidue){
							drawconnections(alphaCIndices[ii], atomInfoNum[ alphaCIndices[ii] ][RESTYN]);
						}
					}
					break;
				case 4:/*connections included with atom types colored*/
    // cerr << "molecule3.4.\n";
							if(!noConnections){
       							FOR(ii,1,lastResidue){
									drawconnections(alphaCIndices[ii], atomInfoNum[ alphaCIndices[ii] ][RESTYN]);
								}
							}
							drawBackbone();break;
				case LASTBBONE:
    // cerr << "molecule3.LAST.\n";
              drawBackbone();break;/* with different colored chains */
			}/*end switch(backBone)*/
     // cerr << "molecule end.\n";
	}


	if(sceneType==KLEIN4){endLune=2;}
	if(sceneType==PARABOLIC||sceneType==MOLECULE){endLune=0;}
  	loxAxisColor;
	FOR(ww,begLune,endLune){/********* The Lune ***********/
		glBegin(GL_LINE_STRIP);
		FOR(jj,0,NAP){ /*********** Axis *****************/
			glVertex3fv(frameAxes[ww][jj]);
		}
		glEnd();
		if(sceneType==LOXODROMIC){
			bandI=4*(int)(loxLength-1)/40 - 2;
		}
		else {bandI=NA-1;}
		FOR(ii,0,bandI){
			if(ii-(int)((float)ii/2)*2==0 && ii<bandI-1){
				glBegin(GL_TRIANGLE_STRIP);
				FOR(jj,0,NAP-1){
				/*luneColor;*/
				switch(sceneType){
						case LOXODROMIC:
							glColor4f(1.0,0.2,1.0,.5);
							break;
						case ELLIPTIC:
							glColor4f(0.2,.25*(float)jj/(NAP-1),.9/*(1-(float)ii/5)*/,0.7);
							break;
						case KLEIN4:
							if(ww==0){glColor4f(0,(1-(float)jj/NAP),(float)ii/(NA-1),0.7);}
							else {glColor4f((float)ii/NA,0.0,(float)jj/NAP,0.7);}
							break;
				}
					v3f(frameLune[ww][0][ii][jj]);
					v3f(frameLune[ww][0][ii+1][jj]);
				}
				glEnd();
			}
		}
	}/* End Lunes*/
	if(sceneType== ELLIPTIC){
			FOR(kk,0,NR){/**** Rings *****/
				FOR(ii,0,5-1){
					if(ii-(int)((float)ii/2)*2==0 && ii<5){
						bgntmesh();
						FOR(jj,0,NAP-1){
							glColor4f(0.0,.5*(float)jj/(NAP-1),(1-(float)kk/5),0.5);
							glVertex3fv(frameRing[0][kk][ii][jj]);
							glVertex3fv(frameRing[0][kk][ii+1][jj]);
						}
						endtmesh();
					}
				}
			}
	}

  // cerr << "done draw.\n";
  dawn(1); 
}
void drawall2(void){/*AKG*/
	int ii, jj, kk, ww/*hichLune*/, begLune=0, endLune=1, bandI;

	/****artNscience: calling the draw routine**************/
	/****drawTheCircle();*********/
if(axx) drawAxes(); if(changeCent) drawSphere(lastCent);/*mark the in-front-of-you spot  */


	FOR(ww,begLune,endLune){/********* The Lune ***********/
		glBegin(GL_LINE_STRIP);
		FOR(jj,0,NAP){ /*********** Axis *****************/
			glVertex3fv(frameAxes[ww][jj]);
		}
		glEnd();

		bandI=4*(int)(loxLength-1)/40 - 2;
		FOR(ii,0,bandI){
			if(ii-(int)((float)ii/2)*2==0 && ii<bandI-1){
				glBegin(GL_TRIANGLE_STRIP);
				FOR(jj,0,NAP-1){
				/*luneColor;*/
				switch(sceneType){
						case LOXODROMIC:
							glColor4f(1.0,0.2,1.0,.5);
							break;
						case ELLIPTIC:
							glColor4f(0.2,.25*(float)jj/(NAP-1),.9/*(1-(float)ii/5)*/,0.7);
							break;
						case KLEIN4:
							if(ww==0){glColor4f(0,(1-(float)jj/NAP),(float)ii/(NA-1),0.7);}
							else {glColor4f((float)ii/NA,0.0,(float)jj/NAP,0.7);}
							break;
				}glColor4f(1.0,0.2,1.0,.5);
					v3f(frameLune[ww][0][ii][jj]);
					v3f(frameLune[ww][0][ii+1][jj]);
				}
				glEnd();
			}
		}
	}/* End Lunes*/

	glLineWidth(3); /*thick*/
			loxPlaneColor;
			glBegin(GL_LINES);
			FOR(ii,0,loxLength){
				glVertex3fv(frameLoxCurvePoints[0][ii]);
			}
			glEnd();
			glBegin(GL_LINES);
			FOR(ii,0,loxLength){
				glVertex3fv(frameLoxCurvePoints[1][ii]);
			}
			glEnd();

}
/**end drawall2**/
/*****Animation********************/
void loxTheScene(int lP, float cF){
	int i, j, a, p, ttt;
	float temp[4][3], temp2[4][3], mult;
	order=6;twist=2*Pi/order;
	FOR(i,0,4){GETS(temp[i],origin);GETS(temp2[i],origin);}

	ttt=(int)lP/10-(int)cF;
	temp2[AAAA][X]=pow(multiplierK,cF)*cos(twist*(cF+ttt));
	temp2[AAAA][Y]=pow(multiplierK,cF)*sin(twist*(cF+ttt));
	temp2[DDDD][X]=1;
	multMatrix(temp,infToF1[0],temp2);

	FOR(a,0,NA){
        	FOR(p,0,NAP){
                	pE(frameLune[0][0][a][p],lune[lastLune][0][a][p],temp);
                }
	}
}
void rotTheScene(float cF){
	int i, j, k, a, p, r, c;
	float temp[4][3], temp2[4][3], twis=2*Pi;
	FOR(i,0,4){GETS(temp[i],origin);GETS(temp2[i],origin);}
	cf[rot][AAAA][X]=cos(twist*cF);
	cf[rot][AAAA][Y]=sin(twist*cF);
	cf[rot][DDDD][X]=1;
	if(!klein2){
		multMatrix(temp,cf[rot],f1ToInf[0]);
		multMatrix(temp2,infToF1[0],temp);
	}
	else {
		multMatrix(temp,cf[rot],f1ToInf[1]);
		multMatrix(temp2,infToF1[1],temp);
	}
	FOR(a,0,NA){
		FOR(p,0,NAP){
			pE(frameLune[0][0][a][p],lune[0][0][a][p],temp2);
		}
	}
	FOR(r,0,NR){
		FOR(c,0,5){
	       		FOR(p,0,NAP){
	       			pE(frameRing[0][r][c][p],
		    		ring[0][0][r][c][p],temp2);
			}
		}
   	}
	switch(pla){
		case 0:
		break;

		case 1:
			FOR(i,0,NH){ /*** Hyperbolic Pencil of Circles ******/
				FOR(j,0,NHAP){
					pE(frameHypPencil[i][j],hypPencil[i][j],temp2);
				}
			}
		break;

		case 2:
			/**/FOR(i,0,NLATS){
	        		FOR(j,0,NAP){
					pE(fPLats[i][j],pLats[i][j],temp2);
				}
			}/**/
			FOR(i,0,NLONGS){
			        FOR(j,0,NAP){
					pE(fPLongs[i][j],pLongs[i][j],temp2);
				}
			}

		break;
	}
	FOR(i,0,NAP){
		GETS(frameAxes[0][p],axes[0][p]);
	}
/**/
	if(sceneType==KLEIN4){
             	FOR(a,0,NA){
                       FOR(p,0,NAP){
                      	   	pE(frameLune[1][0][a][p],
				lune[1][0][a][p],temp2);
                        }
              	}
		FOR(p,0,NAP){
                   	pE(frameAxes[0][p],axes[0][p],temp2);
                   	pE(frameAxes[1][p],axes[1][p],temp2);
                }

	}
}
void movie(){   /* each scene/case makes use of "slides" named still[]... */
	int l/*ayers*/, a/*latitudes*/, p/*oints in each arc*/,
		c/*ircle*/, begTet, endTet, tCount, mC, dd, ddd;
   switch(sceneType){

	case LOXODROMIC:
		if(movieCounter<2*lMax && loxLength<10*NA){
			movieCounter++;loxLength++;loxLengthPerm++;
			loxTheScene(loxLength,(float)movieCounter/lMax);
		}
		else {movieCounter=0;}
	break;
	case ELLIPTIC:
		if((movieCounter/cMax-(float)movieCounter/cMax)!=0){
	 		movieCounter++;
			rotTheScene((float)movieCounter/cMax);
		}
		else {	/* here is a pause */
			if(pauseCounter<pauseMax){pauseCounter++;}
			else {movieCounter++;pauseCounter=0;}
		}
	break;
	case KLEIN4:
		if(movieCounter!=(cMax/2+1) && movieCounter!=cMax+1){
			rotTheScene((float)movieCounter/cMax);
 			movieCounter++;
		}
		else {	/* here is a pause */
			if(pauseCounter<pauseMax) pauseCounter++;
			else {
				pauseCounter=0;
				if(movieCounter!=cMax+1){
					rotTheScene((float)movieCounter/cMax);
				}
				movieCounter++;
			}
		}

		if(movieCounter==cMax+2){
			klein2=(klein2+1)%(2);movieCounter=0;
		}
	break;
	case MOLECULE:
			drawSingleFrame(chainCounter,atomCounter,jointCounter);
			if(jointCounter<spineKnots){
				jointCounter++;
			}
			else {
				jointCounter=0;
				if(atomCounter<lastBB[chainCounter]){
					atomCounter++;
				}
				else {
					atomCounter=0;
					if(chainCounter<lastChain) chainCounter++;
					else chainCounter=1;/*n.b. chains start counting from 1, not 0*/
				}
			}
	break;

   }/* END switch(scen... */


}
/***ANIMEnd*****************************************************************/
/***************************************************************************//**********************************************************************/
void setAtomInfoNumRESTYN(void){/*only needed for initial molecule model*/
	int ii, jj; alphaCIndices[21]=306;
	FOR(ii,1,lastResidue+1+1){ /*residueType[ii]=SER; */

		jj= alphaCIndices[ii]-1;
		while(jj<alphaCIndices[ii+1]-1){
			atomInfoNum[jj][RESTYN]= residueType[ii];
			jj++;
		}

	}/*end FOR loop the residues*/


}
void setAtomRadii(void){
	atomRadius[MARK]=SF*.5;
	atomRadius[HYDR]=SF*.5*(25.0/70.0);
	atomRadius[CARB]=SF*.5;
	atomRadius[NITR]=SF*.5*(65.0/70.0);
	atomRadius[OXYG]=SF*.5*(60.0/70.0);
	atomRadius[SULP]=SF*.5*(60.0/70.0);
}
/* deFault - Sets the default values for state variables.             */
/**********************************************************************/
void reSet(void){
	sceneNum=0;
	movieCounter=0;
	multiplierK=1/.9;
	initializeFigures();
		if(sceneType==MOLECULE){speed=fastSpeed;mov=0;} else{speed=slowSpeed;mov=1;}
	switch(sceneType){
	    case LOXODROMIC:     band=1/**0/**/;sceneLimit=1;loxPath=2;pla=0;
			order=6;twist=2*Pi/order;cMax=30;mov=1;
			loxLength=40;movieCounter=0;
		aff[4*0+0]=0.789819;aff[4*0+1]=0.382883;aff[4*0+2]=-0.479174;aff[4*0+3]=0;
		aff[4*1+0]=-0.599567;aff[4*1+1]=0.31723;aff[4*1+2]=-0.734781;aff[4*1+3]=0;
		aff[4*2+0]=-0.129327;aff[4*2+1]=0.867631;aff[4*2+2]=0.480112;aff[4*2+3]=0;
		aff[4*3+0]=0.0298866;aff[4*3+1]=-0.416265;aff[4*3+2]=-3.42466;aff[4*3+3]=1;
		break;
		case ELLIPTIC: order=5;twist=2*Pi/order;band=1;mov=1;horo=NR;klein2=0;pla=2;
			sceneLimit=order;sceneNum=-1;cMax=30;

		aff[4*0+0]=-0.635972;aff[4*0+1]=-0.143446;aff[4*0+2]=-0.758343;aff[4*0+3]=0;
		aff[4*1+0]=-0.737894;aff[4*1+1]=-0.174933;aff[4*1+2]=0.651926;aff[4*1+3]=0;
		aff[4*2+0]=-0.226161;aff[4*2+1]=0.974124;aff[4*2+2]=0.00540696;aff[4*2+3]=0;
		aff[4*3+0]=0.27603;aff[4*3+1]=-1.15375;aff[4*3+2]=-2.75844;aff[4*3+3]=1;

	    break;
	    case KLEIN4:sceneLimit=8;sceneNum=-1;cMax=64*4;band=1;mov=1;twist=2*Pi;pla=2;
		aff[4*0+0]=0.630936;aff[4*0+1]=-0.0999227;aff[4*0+2]=0.769436;aff[4*0+3]=0;
		aff[4*1+0]=0.765805;aff[4*1+1]=-0.0791756;aff[4*1+2]=-0.638241;aff[4*1+3]=0;
		aff[4*2+0]=0.124688;aff[4*2+1]=0.991876;aff[4*2+2]=0.0265622;aff[4*2+3]=0;
		aff[4*3+0]=0.0298866;aff[4*3+1]=-0.416265;aff[4*3+2]=-3.42466;aff[4*3+3]=1;
 	    break;
		case PARABOLIC:
			pla=2;band=1;
		aff[4*0+0]=0.630936;aff[4*0+1]=-0.0999227;aff[4*0+2]=0.769436;aff[4*0+3]=0;
		aff[4*1+0]=0.765805;aff[4*1+1]=-0.0791756;aff[4*1+2]=-0.638241;aff[4*1+3]=0;
		aff[4*2+0]=0.124688;aff[4*2+1]=0.991876;aff[4*2+2]=0.0265622;aff[4*2+3]=0;
		aff[4*3+0]=0.0298866;aff[4*3+1]=-0.416265;aff[4*3+2]=-3.42466;aff[4*3+3]=1;
 	    break;
 	    case MOLECULE:
			pla=0;
		/**/aff[4*0+0]=0.630936;aff[4*0+1]=-0.0999227;aff[4*0+2]=0.769436;aff[4*0+3]=0;
		aff[4*1+0]=0.765805;aff[4*1+1]=-0.0791756;aff[4*1+2]=-0.638241;aff[4*1+3]=0;
		aff[4*2+0]=0.124688;aff[4*2+1]=0.991876;aff[4*2+2]=0.0265622;aff[4*2+3]=0;
		aff[4*3+0]=0;aff[4*3+1]=0.;aff[4*3+2]=-25.0;aff[4*3+3]=1;
 	    /************************/
		/**/aff2[4*0+0]=0.630936;aff2[4*0+1]=-0.0999227;aff2[4*0+2]=0.769436;aff2[4*0+3]=0;
		aff2[4*1+0]=0.765805;aff2[4*1+1]=-0.0791756;aff2[4*1+2]=-0.638241;aff2[4*1+3]=0;
		aff2[4*2+0]=0.124688;aff2[4*2+1]=0.991876;aff2[4*2+2]=0.0265622;aff2[4*2+3]=0;
		aff2[4*3+0]=7;aff2[4*3+1]=7.;aff2[4*3+2]=-25.0;aff2[4*3+3]=1;
 	    /************************/
 	    break;
	}
		fac1=25;

/******** SZG ***************/
gesture=0; 
eps =0.002; //hand rotation gesture attenuation
del =0.02; //hand translation gesture attenuation
resethand=0;
autotymer(1); /* reset autotymer to start at the beginning */
}
void test100(void){
	int ii, jj;
	strcpy(dummyChars[0][0], "ATOM");
	strcpy(dummyChars[0][1], "ATOM  ");
	dummyNums[0][0]=strcmp(dummyChars[0][0], dummyChars[0][1]);
}
void deFault(void){
th0=5; th1=355;  ta0=5; ta1=355; gap = gap0;
msg=1; nose=.06; mode=TURNMODE;
gesture = 0;
sizz = 1.5;
friz = 0;
if(sceneType==MOLECULE){speed=fastSpeed;} else{speed=slowSpeed;}
torq=.001; focal = 2.; wfar =250.; mysiz=.01; morph=0;mov=0;
aa=1; bb=2; cc=3; dtt=caveyes?DTCAVE:DTCONS;
FOR(ii,0,16) starmat[ii]=aff[ii] =aff[ii] = (ii/4==ii%4);   /* identities */
FOR(ii,0,3)lux[ii]/=NRM(lux); amb = .3; pwr = 10. ;
aff[12]=0; aff[13]= 5; aff[14]= -4.0;   /* place where we can see it */
aff2[12]=0; aff2[13]= 5; aff2[14]= -4.0;   /* place where we can see it */
FOR(ii,0,16) const_aff[ii] = const_aff2[ii]= (ii/4==ii%4);   /* identity */
const_aff[13]= 5.0;   const_aff2[14]= -5.0;/* place where we can see it */

autotymer(1); /* reset autotymer to start at the beginning */

tt=0;
cockpit=0;
cockpit = 0;
autopilot=0;
showslate = 0;


speedMed=.01;
geo=0;
/*akgeowall*/
	if(geo){binoc=-1;}/*n.b. if(binoc) returns "true" for any binoc value other than 0*/

	else{binoc=0;} /*begin with binocular vision turned off*/

	mov=0;lMax=10;pauseMax=1;
	multiplierK=1/.9;
	order=6;twist=2*Pi/order;lastLune=9;
	fixedPnts[0][0][X]=1;fixedPnts[0][0][Y]=0;fixedPnts[0][0][Z]=0;
	fixedPnts[0][1][X]=-1;fixedPnts[0][1][Y]=0;fixedPnts[0][1][Z]=0;
	fixedPnts[1][0][X]=0;fixedPnts[1][0][Y]=1;fixedPnts[1][0][Z]=0;
	fixedPnts[1][1][X]=0;fixedPnts[1][1][Y]=-1;fixedPnts[1][1][Z]=0;
	baseVertices[1][X]=1;baseVertices[2][X]=.5;baseVertices[2][Y]=pow(3.f,.5f)/2;

	FOR(ii,0,TH){
		theta[ii]=Pi/6*(1-(float)ii/TH*.8);
	}

/**artNscience: colors***********/
col[101][0]=1; col[101][1]=0; col[101][2]=0;
col[102][0]=0; col[102][1]=1; col[102][2]=0;
col[103][0]=1; col[103][1]=1; col[103][2]=0;
col[104][0]=0; col[103][1]=0; col[104][2]=1;
col[105][0]=1; col[105][1]=0; col[105][2]=1;
col[106][0]=0; col[106][1]=1; col[106][2]=1;
col[107][0]=1; col[107][1]=1; col[107][2]=1;
FOR(kk,1,5){
FOR(ii,1,8){
	FOR(jj,0,3){col[107+kk*ii][jj]=col[100+ii][jj];}
}}
	/**/ col[BKBN][0]=0.8;col[BKBN][1]=0.8;col[BKBN][2]=0.8;	/*backbone is lt gray*/
	/**/ col[MARK][0]=1.0;col[MARK][1]=0.2;col[MARK][2]=1.0;	/*marker is purple*/
	/**/ col[HYDR][0]=0.7;col[HYDR][1]=0.7;col[HYDR][2]=1.0;	/*hydrogen is blue-white*/
	/**/ col[CARB][0]=0.7;col[CARB][1]=0.7;col[CARB][2]=0.7;	/*carbon is gray*/
	/**/ col[NITR][0]=0.3;col[NITR][1]=0.3;col[NITR][2]=1.0;	/*nitrogen is blue*/
	/**/ col[OXYG][0]=1.0;col[OXYG][1]=0.3;col[OXYG][2]=0.3;	/*oxygen is red*/
	/**/ col[SULP][0]=1.0;col[SULP][1]=1.0;col[SULP][2]=0.3;	/*sulphur is yellow*/
	/**/ col[PHIL][0]=0.5;col[PHIL][1]=1.0;col[PHIL][2]=1.0;	/*hydrophilic is cyan*/
	/**/ col[PHOB][0]=0.7;col[PHOB][1]=0.7;col[PHOB][2]=0.3;	/*hydrophobic is brown*/
	/**/ col[CHRG][0]=0.3;col[CHRG][1]=1.0;col[CHRG][2]=1.0;	/*charged (or all) is green*/
	/**/ col[CHRGP][0]=0.5;col[CHRGP][1]=0.5;col[CHRGP][2]=1.0;	/*pos charged is lt blue*/
	/**/ col[CHRGN][0]=1.0;col[CHRGN][1]=0.3;col[CHRGN][2]=0.3;	/*neg charged is lt red*/
	/**/ col[ELSE][0]=0.8;col[ELSE][1]=0.8;col[ELSE][2]=0.8;	/*else is lt gray*/

/**/aCent[0][X]=-8.901*SF; aCent[0][Y]=4.127*SF;  aCent[0][Z]=-0.555*SF; aCent[0][W]=0;
aCent[1][X]=-8.901*SF; aCent[1][Y]=4.127*SF;  aCent[1][Z]=-0.555*SF; aCent[1][W]=NITR;
aCent[2][X]=-8.608*SF; aCent[2][Y]=3.135*SF;  aCent[2][Z]=-1.618*SF; aCent[2][W]=CARB;
aCent[3][X]=-7.117*SF; aCent[3][Y]=2.964*SF;  aCent[3][Z]=-1.897*SF; aCent[3][W]=CARB;
aCent[4][X]=-6.634*SF; aCent[4][Y]=1.849*SF;  aCent[4][Z]=-1.758*SF; aCent[4][W]=OXYG;
aCent[5][X]=-9.437*SF; aCent[5][Y]=3.396*SF;  aCent[5][Z]=-2.889*SF; aCent[5][W]=CARB;
aCent[6][X]=-10.915*SF; aCent[6][Y]=3.130*SF;  aCent[6][Z]=-2.611*SF; aCent[6][W]=CARB;
aCent[7][X]=-11.269*SF; aCent[7][Y]=2.700*SF;  aCent[7][Z]=-1.524*SF; aCent[7][W]=OXYG;
aCent[8][X]=-11.806*SF; aCent[8][Y]=3.406*SF;  aCent[8][Z]=-3.543*SF; aCent[8][W]=NITR;
aCent[9][X]=-8.330*SF; aCent[9][Y]=3.957*SF;  aCent[9][Z]=0.261*SF; aCent[9][W]=HYDR;
aCent[10][X]=-8.740*SF; aCent[10][Y]=5.068*SF;  aCent[10][Z]=-0.889*SF; aCent[10][W]=HYDR;
aCent[11][X]=-9.877*SF; aCent[11][Y]=4.041*SF;  aCent[11][Z]=-0.293*SF; aCent[11][W]=HYDR;
aCent[12][X]=-8.930*SF; aCent[12][Y]=2.162*SF;  aCent[12][Z]=-1.239*SF; aCent[12][W]=HYDR;
aCent[13][X]=-9.310*SF; aCent[13][Y]=4.417*SF;  aCent[13][Z]=-3.193*SF; aCent[13][W]=HYDR;
aCent[14][X]=-9.108*SF; aCent[14][Y]=2.719*SF;  aCent[14][Z]=-3.679*SF; aCent[14][W]=HYDR;
aCent[15][X]=-11.572*SF; aCent[15][Y]=3.791*SF;  aCent[15][Z]=-4.444*SF; aCent[15][W]=HYDR;
aCent[16][X]=-12.757*SF; aCent[16][Y]=3.183*SF;  aCent[16][Z]=-3.294*SF; aCent[16][W]=HYDR;
aCent[17][X]=-6.379*SF; aCent[17][Y]=4.031*SF;  aCent[17][Z]=-2.228*SF; aCent[17][W]=NITR;
aCent[18][X]=-4.923*SF; aCent[18][Y]=4.002*SF;  aCent[18][Z]=-2.452*SF; aCent[18][W]=CARB;
aCent[19][X]=-4.136*SF; aCent[19][Y]=3.187*SF;  aCent[19][Z]=-1.404*SF; aCent[19][W]=CARB;
aCent[20][X]=-3.391*SF; aCent[20][Y]=2.274*SF; aCent[20][Z]=-1.760*SF; aCent[20][W]=OXYG;
aCent[21][X]=-4.411*SF; aCent[21][Y]=5.450*SF; aCent[21][Z]=-2.619*SF; aCent[21][W]=CARB;
aCent[22][X]=-4.795*SF; aCent[22][Y]=6.450*SF; aCent[22][Z]=-1.495*SF; aCent[22][W]=CARB;
aCent[23][X]=-3.612*SF; aCent[23][Y]=6.803*SF; aCent[23][Z]=-0.599*SF; aCent[23][W]=CARB;
aCent[24][X]=-5.351*SF; aCent[24][Y]=7.748*SF; aCent[24][Z]=-2.084*SF; aCent[24][W]=CARB;
aCent[25][X]=-6.821*SF; aCent[25][Y]=4.923*SF; aCent[25][Z]=-2.394*SF; aCent[25][W]=HYDR;
aCent[26][X]=-4.750*SF; aCent[26][Y]=3.494*SF; aCent[26][Z]=-3.403*SF; aCent[26][W]=HYDR;
aCent[27][X]=-3.340*SF; aCent[27][Y]=5.414*SF; aCent[27][Z]=-2.672*SF; aCent[27][W]=HYDR;
aCent[28][X]=-4.813*SF; aCent[28][Y]=5.817*SF; aCent[28][Z]=-3.564*SF; aCent[28][W]=HYDR;
aCent[29][X]=-5.568*SF; aCent[29][Y]=6.022*SF; aCent[29][Z]=-0.858*SF; aCent[29][W]=HYDR;
aCent[30][X]=-3.207*SF; aCent[30][Y]=5.905*SF; aCent[30][Z]=-0.146*SF; aCent[30][W]=HYDR;
aCent[31][X]=-2.841*SF; aCent[31][Y]=7.304*SF; aCent[31][Z]=-1.183*SF; aCent[31][W]=HYDR;
aCent[32][X]=-3.929*SF; aCent[32][Y]=7.477*SF; aCent[32][Z]=0.197*SF; aCent[32][W]=HYDR;
aCent[33][X]=-4.607*SF; aCent[33][Y]=8.209*SF; aCent[33][Z]=-2.736*SF; aCent[33][W]=HYDR;
aCent[34][X]=-6.255*SF; aCent[34][Y]=7.544*SF; aCent[34][Z]=-2.657*SF; aCent[34][W]=HYDR;
aCent[35][X]=-5.592*SF; aCent[35][Y]=8.445*SF; aCent[35][Z]=-1.281*SF; aCent[35][W]=HYDR;
aCent[36][X]=-4.354*SF; aCent[36][Y]=3.455*SF; aCent[36][Z]=-0.111*SF; aCent[36][W]=NITR;
aCent[37][X]=-3.690*SF; aCent[37][Y]=2.738*SF; aCent[37][Z]=0.981*SF; aCent[37][W]=CARB;
aCent[38][X]=-4.102*SF; aCent[38][Y]=1.256*SF; aCent[38][Z]=1.074*SF; aCent[38][W]=CARB;
aCent[39][X]=-3.291*SF; aCent[39][Y]=0.409*SF; aCent[39][Z]=1.442*SF; aCent[39][W]=OXYG;
aCent[40][X]=-3.964*SF; aCent[40][Y]=3.472*SF; aCent[40][Z]=2.302*SF; aCent[40][W]=CARB;
aCent[41][X]=-2.824*SF; aCent[41][Y]=3.339*SF; aCent[41][Z]=3.290*SF; aCent[41][W]=CARB;
aCent[42][X]=-1.820*SF; aCent[42][Y]=4.326*SF; aCent[42][Z]=3.332*SF; aCent[42][W]=CARB;
aCent[43][X]=-2.746*SF; aCent[43][Y]=2.217*SF; aCent[43][Z]=4.138*SF; aCent[43][W]=CARB;
aCent[44][X]=-0.725*SF; aCent[44][Y]=4.185*SF; aCent[44][Z]=4.205*SF; aCent[44][W]=CARB;
aCent[45][X]=-1.657*SF; aCent[45][Y]=2.076*SF; aCent[45][Z]=5.018*SF; aCent[45][W]=CARB;
aCent[46][X]=-0.639*SF; aCent[46][Y]=3.053*SF; aCent[46][Z]=5.043*SF; aCent[46][W]=CARB;
aCent[47][X]=0.433*SF; aCent[47][Y]=2.881*SF; aCent[47][Z]=5.861*SF; aCent[47][W]=OXYG;
aCent[48][X]=-4.934*SF; aCent[48][Y]=4.245*SF; aCent[48][Z]=0.120*SF; aCent[48][W]=HYDR;
aCent[49][X]=-2.615*SF; aCent[49][Y]=2.768*SF; aCent[49][Z]=0.796*SF; aCent[49][W]=HYDR;
aCent[50][X]=-4.117*SF; aCent[50][Y]=4.513*SF; aCent[50][Z]=2.091*SF; aCent[50][W]=HYDR;
aCent[51][X]=-4.886*SF; aCent[51][Y]=3.096*SF; aCent[51][Z]=2.750*SF; aCent[51][W]=HYDR;
aCent[52][X]=-1.877*SF; aCent[52][Y]=5.200*SF; aCent[52][Z]=2.695*SF; aCent[52][W]=HYDR;
aCent[53][X]=-3.513*SF; aCent[53][Y]=1.456*SF; aCent[53][Z]=4.101*SF; aCent[53][W]=HYDR;
aCent[54][X]=0.033*SF; aCent[54][Y]=4.952*SF; aCent[54][Z]=4.233*SF; aCent[54][W]=HYDR;
aCent[55][X]=-1.576*SF; aCent[55][Y]=1.221*SF; aCent[55][Z]=5.669*SF; aCent[55][W]=HYDR;
aCent[56][X]=1.187*SF; aCent[56][Y]=3.395*SF; aCent[56][Z]=5.567*SF; aCent[56][W]=HYDR;
aCent[57][X]=-5.342*SF; aCent[57][Y]=0.925*SF; aCent[57][Z]=0.689*SF; aCent[57][W]=NITR;
aCent[58][X]=-5.857*SF; aCent[58][Y]=-0.449*SF; aCent[58][Z]=0.613*SF; aCent[58][W]=CARB;
aCent[59][X]=-5.089*SF; aCent[59][Y]=-1.221*SF; aCent[59][Z]=-0.470*SF; aCent[59][W]=CARB;
aCent[60][X]=-4.621*SF; aCent[60][Y]=-2.334*SF; aCent[60][Z]=-0.226*SF; aCent[60][W]=OXYG;
aCent[61][X]=-7.386*SF; aCent[61][Y]=-0.466*SF; aCent[61][Z]=0.343*SF; aCent[61][W]=CARB;
aCent[62][X]=-8.197*SF; aCent[62][Y]=0.540*SF; aCent[62][Z]=1.197*SF; aCent[62][W]=CARB;
aCent[63][X]=-7.959*SF; aCent[63][Y]=-1.884*SF; aCent[63][Z]=0.501*SF; aCent[63][W]=CARB;
aCent[64][X]=-8.019*SF; aCent[64][Y]=0.412*SF; aCent[64][Z]=2.715*SF; aCent[64][W]=CARB;
aCent[65][X]=-5.906*SF; aCent[65][Y]=1.656*SF; aCent[65][Z]=0.283*SF; aCent[65][W]=HYDR;
aCent[66][X]=-5.670*SF; aCent[66][Y]=-0.941*SF; aCent[66][Z]=1.568*SF; aCent[66][W]=HYDR;
aCent[67][X]=-7.554*SF; aCent[67][Y]=-0.192*SF; aCent[67][Z]=-0.697*SF; aCent[67][W]=HYDR;
aCent[68][X]=-7.900*SF; aCent[68][Y]=1.531*SF; aCent[68][Z]=0.912*SF; aCent[68][W]=HYDR;
aCent[69][X]=-9.257*SF; aCent[69][Y]=0.424*SF; aCent[69][Z]=0.964*SF; aCent[69][W]=HYDR;
aCent[70][X]=-7.509*SF; aCent[70][Y]=-2.555*SF; aCent[70][Z]=-0.232*SF; aCent[70][W]=HYDR;
aCent[71][X]=-7.759*SF; aCent[71][Y]=-2.271*SF; aCent[71][Z]=1.501*SF; aCent[71][W]=HYDR;
aCent[72][X]=-9.036*SF; aCent[72][Y]=-1.871*SF; aCent[72][Z]=0.332*SF; aCent[72][W]=HYDR;
aCent[73][X]=-8.306*SF; aCent[73][Y]=-0.585*SF; aCent[73][Z]=3.049*SF; aCent[73][W]=HYDR;
aCent[74][X]=-6.983*SF; aCent[74][Y]=0.606*SF; aCent[74][Z]=2.995*SF; aCent[74][W]=HYDR;
aCent[75][X]=-8.656*SF; aCent[75][Y]=1.144*SF; aCent[75][Z]=3.213*SF; aCent[75][W]=HYDR;
aCent[76][X]=-4.907*SF; aCent[76][Y]=-0.601*SF; aCent[76][Z]=-1.645*SF; aCent[76][W]=NITR;
aCent[77][X]=-4.122*SF; aCent[77][Y]=-1.167*SF; aCent[77][Z]=-2.743*SF; aCent[77][W]=CARB;
aCent[78][X]=-2.629*SF; aCent[78][Y]=-1.321*SF; aCent[78][Z]=-2.390*SF; aCent[78][W]=CARB;
aCent[79][X]=-1.986*SF; aCent[79][Y]=-2.240*SF; aCent[79][Z]=-2.884*SF; aCent[79][W]=OXYG;
aCent[80][X]=-4.292*SF; aCent[80][Y]=-0.313*SF; aCent[80][Z]=-4.013*SF; aCent[80][W]=CARB;
aCent[81][X]=-4.244*SF; aCent[81][Y]=-1.171*SF; aCent[81][Z]=-5.290*SF; aCent[81][W]=CARB;
aCent[82][X]=-5.576*SF; aCent[82][Y]=-1.860*SF; aCent[82][Z]=-5.585*SF; aCent[82][W]=CARB;
aCent[83][X]=-5.769*SF; aCent[83][Y]=-3.044*SF; aCent[83][Z]=-5.335*SF; aCent[83][W]=OXYG;
aCent[84][X]=-6.532*SF; aCent[84][Y]=-1.146*SF; aCent[84][Z]=-6.152*SF; aCent[84][W]=NITR;
aCent[85][X]=-5.327*SF; aCent[85][Y]=0.318*SF; aCent[85][Z]=-1.763*SF; aCent[85][W]=HYDR;
aCent[86][X]=-4.517*SF; aCent[86][Y]=-2.162*SF; aCent[86][Z]=-2.940*SF; aCent[86][W]=HYDR;
aCent[87][X]=-5.238*SF; aCent[87][Y]=0.191*SF; aCent[87][Z]=-3.969*SF; aCent[87][W]=HYDR;
aCent[88][X]=-3.492*SF; aCent[88][Y]=0.429*SF; aCent[88][Z]=-4.053*SF; aCent[88][W]=HYDR;
aCent[89][X]=-3.993*SF; aCent[89][Y]=-0.539*SF; aCent[89][Z]=-6.120*SF; aCent[89][W]=HYDR;
aCent[90][X]=-3.458*SF; aCent[90][Y]=-1.923*SF; aCent[90][Z]=-5.205*SF; aCent[90][W]=HYDR;
aCent[91][X]=-6.389*SF; aCent[91][Y]=-0.184*SF; aCent[91][Z]=-6.408*SF; aCent[91][W]=HYDR;
aCent[92][X]=-7.392*SF; aCent[92][Y]=-1.635*SF; aCent[92][Z]=-6.335*SF; aCent[92][W]=HYDR;
aCent[93][X]=-2.074*SF; aCent[93][Y]=-0.459*SF; aCent[93][Z]=-1.528*SF; aCent[93][W]=NITR;
aCent[94][X]=-0.716*SF; aCent[94][Y]=-0.631*SF; aCent[94][Z]=-0.993*SF; aCent[94][W]=CARB;
aCent[95][X]=-0.631*SF; aCent[95][Y]=-1.766*SF; aCent[95][Z]=0.044*SF; aCent[95][W]=CARB;
aCent[96][X]=0.295*SF; aCent[96][Y]=-2.579*SF; aCent[96][Z]=-0.004*SF; aCent[96][W]=OXYG;
aCent[97][X]=-0.221*SF; aCent[97][Y]=0.703*SF; aCent[97][Z]=-0.417*SF; aCent[97][W]=CARB;
aCent[98][X]=1.148*SF; aCent[98][Y]=0.652*SF; aCent[98][Z]=0.194*SF; aCent[98][W]=CARB;
aCent[99][X]=2.319*SF; aCent[99][Y]=0.664*SF; aCent[99][Z]=-0.482*SF; aCent[99][W]=CARB;
aCent[100][X]=1.508*SF;aCent[100][Y]=0.564*SF;aCent[100][Z]=1.606*SF;aCent[100][W]=CARB;
aCent[101][X]=3.371*SF;aCent[101][Y]=0.560*SF;aCent[101][Z]=0.411*SF;aCent[101][W]=NITR;
aCent[102][X]=2.928*SF;aCent[102][Y]=0.515*SF;aCent[102][Z]=1.710*SF;aCent[102][W]=CARB;
aCent[103][X]=0.779*SF;aCent[103][Y]=0.524*SF;aCent[103][Z]=2.812*SF;aCent[103][W]=CARB;
aCent[104][X]=3.599*SF;aCent[104][Y]=0.445*SF;aCent[104][Z]=2.938*SF;aCent[104][W]=CARB;
aCent[105][X]=1.439*SF;aCent[105][Y]=0.433*SF;aCent[105][Z]=4.053*SF;aCent[105][W]=CARB;
aCent[106][X]=2.842*SF;aCent[106][Y]=0.407*SF;aCent[106][Z]=4.120*SF;aCent[106][W]=CARB;
aCent[107][X]=-2.624*SF;aCent[107][Y]=0.343*SF;aCent[107][Z]=-1.242*SF;aCent[107][W]=HYDR;
aCent[108][X]=-0.052*SF;aCent[108][Y]=-0.908*SF;aCent[108][Z]=-1.813*SF;aCent[108][W]=HYDR;
aCent[109][X]=-0.206*SF;aCent[109][Y]=1.425*SF;aCent[109][Z]=-1.211*SF;aCent[109][W]=HYDR;
aCent[110][X]=-0.921*SF;aCent[110][Y]=1.044*SF;aCent[110][Z]=0.344*SF;aCent[110][W]=HYDR;
aCent[111][X]=2.412*SF;aCent[111][Y]=0.733*SF;aCent[111][Z]=-1.558*SF;aCent[111][W]=HYDR;
aCent[112][X]=4.360*SF;aCent[112][Y]=0.536*SF;aCent[112][Z]=0.156*SF;aCent[112][W]=HYDR;
aCent[113][X]=-0.299*SF;aCent[113][Y]=0.571*SF;aCent[113][Z]=2.773*SF;aCent[113][W]=HYDR;
aCent[114][X]=4.679*SF;aCent[114][Y]=0.418*SF;aCent[114][Z]=2.961*SF;aCent[114][W]=HYDR;
aCent[115][X]=0.862*SF;aCent[115][Y]=0.400*SF;aCent[115][Z]=4.966*SF;aCent[115][W]=HYDR;
aCent[116][X]=3.334*SF;aCent[116][Y]=0.360*SF;aCent[116][Z]=5.081*SF;aCent[116][W]=HYDR;
aCent[117][X]=-1.600*SF;aCent[117][Y]=-1.860*SF;aCent[117][Z]=0.967*SF;aCent[117][W]=NITR;
aCent[118][X]=-1.641*SF;aCent[118][Y]=-2.932*SF;aCent[118][Z]=1.963*SF;aCent[118][W]=CARB;
aCent[119][X]=-1.847*SF;aCent[119][Y]=-4.319*SF;aCent[119][Z]=1.342*SF;aCent[119][W]=CARB;
aCent[120][X]=-1.144*SF;aCent[120][Y]=-5.248*SF;aCent[120][Z]=1.742*SF;aCent[120][W]=OXYG;
aCent[121][X]=-2.710*SF;aCent[121][Y]=-2.645*SF;aCent[121][Z]=3.033*SF;aCent[121][W]=CARB;
aCent[122][X]=-2.301*SF;aCent[122][Y]=-1.579*SF;aCent[122][Z]=4.069*SF;aCent[122][W]=CARB;
aCent[123][X]=-3.475*SF;aCent[123][Y]=-1.323*SF;aCent[123][Z]=5.018*SF;aCent[123][W]=CARB;
aCent[124][X]=-1.093*SF;aCent[124][Y]=-2.007*SF;aCent[124][Z]=4.914*SF;aCent[124][W]=CARB;
aCent[125][X]=-2.316*SF;aCent[125][Y]=-1.137*SF;aCent[125][Z]=0.994*SF;aCent[125][W]=HYDR;
aCent[126][X]=-0.666*SF;aCent[126][Y]=-2.978*SF;aCent[126][Z]=2.445*SF;aCent[126][W]=HYDR;
aCent[127][X]=-3.600*SF;aCent[127][Y]=-2.308*SF;aCent[127][Z]=2.537*SF;aCent[127][W]=HYDR;
aCent[128][X]=-2.921*SF;aCent[128][Y]=-3.571*SF;aCent[128][Z]=3.572*SF;aCent[128][W]=HYDR;
aCent[129][X]=-2.061*SF;aCent[129][Y]=-0.649*SF;aCent[129][Z]=3.560*SF;aCent[129][W]=HYDR;
aCent[130][X]=-4.343*SF;aCent[130][Y]=-0.992*SF;aCent[130][Z]=4.449*SF;aCent[130][W]=HYDR;
aCent[131][X]=-3.725*SF;aCent[131][Y]=-2.237*SF;aCent[131][Z]=5.560*SF;aCent[131][W]=HYDR;
aCent[132][X]=-3.211*SF;aCent[132][Y]=-0.549*SF;aCent[132][Z]=5.739*SF;aCent[132][W]=HYDR;
aCent[133][X]=-1.270*SF;aCent[133][Y]=-2.989*SF;aCent[133][Z]=5.354*SF;aCent[133][W]=HYDR;
aCent[134][X]=-0.195*SF;aCent[134][Y]=-2.045*SF;aCent[134][Z]=4.300*SF;aCent[134][W]=HYDR;
aCent[135][X]=-0.922*SF;aCent[135][Y]=-1.286*SF;aCent[135][Z]=5.712*SF;aCent[135][W]=HYDR;
aCent[136][X]=-2.753*SF;aCent[136][Y]=-4.481*SF;aCent[136][Z]=0.360*SF;aCent[136][W]=NITR;
aCent[137][X]=-3.024*SF;aCent[137][Y]=-5.791*SF;aCent[137][Z]=-0.269*SF;aCent[137][W]=CARB;
aCent[138][X]=-1.796*SF;aCent[138][Y]=-6.427*SF;aCent[138][Z]=-0.937*SF;aCent[138][W]=CARB;
aCent[139][X]=-1.719*SF;aCent[139][Y]=-7.648*SF;aCent[139][Z]=-1.030*SF;aCent[139][W]=OXYG;
aCent[140][X]=-4.224*SF;aCent[140][Y]=-5.697*SF;aCent[140][Z]=-1.232*SF;aCent[140][W]=CARB;
aCent[141][X]=-3.930*SF;aCent[141][Y]=-5.009*SF;aCent[141][Z]=-2.577*SF;aCent[141][W]=CARB;
aCent[142][X]=-3.682*SF;aCent[142][Y]=-5.986*SF;aCent[142][Z]=-3.736*SF;aCent[142][W]=CARB;
aCent[143][X]=-3.494*SF;aCent[143][Y]=-5.199*SF;aCent[143][Z]=-5.039*SF;aCent[143][W]=CARB;
aCent[144][X]=-4.563*SF;aCent[144][Y]=-5.483*SF;aCent[144][Z]=-6.023*SF;aCent[144][W]=NITR;
aCent[145][X]=-3.321*SF;aCent[145][Y]=-3.675*SF;aCent[145][Z]=0.097*SF;aCent[145][W]=HYDR;
aCent[146][X]=-3.309*SF;aCent[146][Y]=-6.478*SF;aCent[146][Z]=0.529*SF;aCent[146][W]=HYDR;
aCent[147][X]=-4.565*SF;aCent[147][Y]=-6.694*SF;aCent[147][Z]=-1.436*SF;aCent[147][W]=HYDR;
aCent[148][X]=-5.019*SF;aCent[148][Y]=-5.143*SF;aCent[148][Z]=-0.731*SF;aCent[148][W]=HYDR;
aCent[149][X]=-4.769*SF;aCent[149][Y]=-4.390*SF;aCent[149][Z]=-2.830*SF;aCent[149][W]=HYDR;
aCent[150][X]=-3.062*SF;aCent[150][Y]=-4.368*SF;aCent[150][Z]=-2.469*SF;aCent[150][W]=HYDR;
aCent[151][X]=-2.799*SF;aCent[151][Y]=-6.562*SF;aCent[151][Z]=-3.536*SF;aCent[151][W]=HYDR;
aCent[152][X]=-4.524*SF;aCent[152][Y]=-6.674*SF;aCent[152][Z]=-3.818*SF;aCent[152][W]=HYDR;
aCent[153][X]=-3.502*SF;aCent[153][Y]=-4.150*SF;aCent[153][Z]=-4.813*SF;aCent[153][W]=HYDR;
aCent[154][X]=-2.511*SF;aCent[154][Y]=-5.439*SF;aCent[154][Z]=-5.457*SF;aCent[154][W]=HYDR;
aCent[155][X]=-4.621*SF;aCent[155][Y]=-6.474*SF;aCent[155][Z]=-6.211*SF;aCent[155][W]=HYDR;
aCent[156][X]=-5.442*SF;aCent[156][Y]=-5.124*SF;aCent[156][Z]=-5.657*SF;aCent[156][W]=HYDR;
aCent[157][X]=-4.382*SF;aCent[157][Y]=-4.983*SF;aCent[157][Z]=-6.881*SF;aCent[157][W]=HYDR;
aCent[158][X]=-0.828*SF;aCent[158][Y]=-5.607*SF;aCent[158][Z]=-1.355*SF;aCent[158][W]=NITR;
aCent[159][X]=0.466*SF;aCent[159][Y]=-6.016*SF;aCent[159][Z]=-1.905*SF;aCent[159][W]=CARB;
aCent[160][X]=1.481*SF;aCent[160][Y]=-6.464*SF;aCent[160][Z]=-0.832*SF;aCent[160][W]=CARB;
aCent[161][X]=2.545*SF;aCent[161][Y]=-6.971*SF;aCent[161][Z]=-1.194*SF;aCent[161][W]=OXYG;
aCent[162][X]=1.033*SF;aCent[162][Y]=-4.839*SF;aCent[162][Z]=-2.724*SF;aCent[162][W]=CARB;
aCent[163][X]=0.672*SF;aCent[163][Y]=-4.906*SF;aCent[163][Z]=-4.210*SF;aCent[163][W]=CARB;
aCent[164][X]=-0.532*SF;aCent[164][Y]=-5.051*SF;aCent[164][Z]=-4.522*SF;aCent[164][W]=OXYG;
aCent[165][X]=1.627*SF;aCent[165][Y]=-4.815*SF;aCent[165][Z]=-5.017*SF;aCent[165][W]=OXYG;
aCent[166][X]=-1.010*SF;aCent[166][Y]=-4.616*SF;aCent[166][Z]=-1.291*SF;aCent[166][W]=HYDR;
aCent[167][X]=0.319*SF;aCent[167][Y]=-6.867*SF;aCent[167][Z]=-2.574*SF;aCent[167][W]=HYDR;
aCent[168][X]=0.644*SF;aCent[168][Y]=-3.924*SF;aCent[168][Z]=-2.320*SF;aCent[168][W]=HYDR;
aCent[169][X]=2.116*SF;aCent[169][Y]=-4.837*SF;aCent[169][Z]=-2.650*SF;aCent[169][W]=HYDR;
aCent[170][X]=1.185*SF;aCent[170][Y]=-6.278*SF;aCent[170][Z]=0.464*SF;aCent[170][W]=NITR;
aCent[171][X]=2.060*SF;aCent[171][Y]=-6.618*SF;aCent[171][Z]=1.593*SF;aCent[171][W]=CARB;
aCent[172][X]=2.628*SF;aCent[172][Y]=-5.412*SF;aCent[172][Z]=2.353*SF;aCent[172][W]=CARB;
aCent[173][X]=3.496*SF;aCent[173][Y]=-5.594*SF;aCent[173][Z]=3.208*SF;aCent[173][W]=OXYG;
aCent[174][X]=0.265*SF;aCent[174][Y]=-5.908*SF;aCent[174][Z]=0.693*SF;aCent[174][W]=HYDR;
aCent[175][X]=1.486*SF;aCent[175][Y]=-7.214*SF;aCent[175][Z]=2.304*SF;aCent[175][W]=HYDR;
aCent[176][X]=2.897*SF;aCent[176][Y]=-7.228*SF;aCent[176][Z]=1.252*SF;aCent[176][W]=HYDR;
aCent[177][X]=2.172*SF;aCent[177][Y]=-4.187*SF;aCent[177][Z]=2.055*SF;aCent[177][W]=NITR;
aCent[178][X]=2.626*SF;aCent[178][Y]=-2.967*SF;aCent[178][Z]=2.723*SF;aCent[178][W]=CARB;
aCent[179][X]=4.157*SF;aCent[179][Y]=-2.802*SF;aCent[179][Z]=2.654*SF;aCent[179][W]=CARB;
aCent[180][X]=4.710*SF;aCent[180][Y]=-2.829*SF;aCent[180][Z]=1.551*SF;aCent[180][W]=OXYG;
aCent[181][X]=1.481*SF;aCent[181][Y]=-4.089*SF;aCent[181][Z]=1.319*SF;aCent[181][W]=HYDR;
aCent[182][X]=2.164*SF;aCent[182][Y]=-2.109*SF;aCent[182][Z]=2.237*SF;aCent[182][W]=HYDR;
aCent[183][X]=2.280*SF;aCent[183][Y]=-2.997*SF;aCent[183][Z]=3.753*SF;aCent[183][W]=HYDR;
aCent[184][X]=4.871*SF;aCent[184][Y]=-2.651*SF;aCent[184][Z]=3.794*SF;aCent[184][W]=NITR;
aCent[185][X]=6.333*SF;aCent[185][Y]=-2.533*SF;aCent[185][Z]=3.806*SF;aCent[185][W]=CARB;
aCent[186][X]=7.058*SF;aCent[186][Y]=-3.729*SF;aCent[186][Z]=3.165*SF;aCent[186][W]=CARB;
aCent[187][X]=8.139*SF;aCent[187][Y]=-3.562*SF;aCent[187][Z]=2.601*SF;aCent[187][W]=OXYG;
aCent[188][X]=6.740*SF;aCent[188][Y]=-2.387*SF;aCent[188][Z]=5.279*SF;aCent[188][W]=CARB;
aCent[189][X]=5.460*SF;aCent[189][Y]=-1.952*SF;aCent[189][Z]=5.987*SF;aCent[189][W]=CARB;
aCent[190][X]=4.362*SF;aCent[190][Y]=-2.615*SF;aCent[190][Z]=5.160*SF;aCent[190][W]=CARB;
aCent[191][X]=6.611*SF;aCent[191][Y]=-1.626*SF;aCent[191][Z]=3.267*SF;aCent[191][W]=HYDR;
aCent[192][X]=7.091*SF;aCent[192][Y]=-3.323*SF;aCent[192][Z]=5.670*SF;aCent[192][W]=HYDR;
aCent[193][X]=7.531*SF;aCent[193][Y]=-1.647*SF;aCent[193][Z]=5.403*SF;aCent[193][W]=HYDR;
aCent[194][X]=5.443*SF;aCent[194][Y]=-2.302*SF;aCent[194][Z]=7.001*SF;aCent[194][W]=HYDR;
aCent[195][X]=5.358*SF;aCent[195][Y]=-0.867*SF;aCent[195][Z]=5.929*SF;aCent[195][W]=HYDR;
aCent[196][X]=4.173*SF;aCent[196][Y]=-3.609*SF;aCent[196][Z]=5.516*SF;aCent[196][W]=HYDR;
aCent[197][X]=3.440*SF;aCent[197][Y]=-2.042*SF;aCent[197][Z]=5.246*SF;aCent[197][W]=HYDR;
aCent[198][X]=6.463*SF;aCent[198][Y]=-4.929*SF;aCent[198][Z]=3.205*SF;aCent[198][W]=NITR;
aCent[199][X]=7.049*SF;aCent[199][Y]=-6.179*SF;aCent[199][Z]=2.704*SF;aCent[199][W]=CARB;
aCent[200][X]=6.897*SF;aCent[200][Y]=-6.369*SF;aCent[200][Z]=1.185*SF;aCent[200][W]=CARB;
aCent[201][X]=7.025*SF;aCent[201][Y]=-7.488*SF;aCent[201][Z]=0.697*SF;aCent[201][W]=OXYG;
aCent[202][X]=6.458*SF;aCent[202][Y]=-7.371*SF;aCent[202][Z]=3.472*SF;aCent[202][W]=CARB;
aCent[203][X]=6.763*SF;aCent[203][Y]=-7.264*SF;aCent[203][Z]=4.850*SF;aCent[203][W]=OXYG;
aCent[204][X]=5.535*SF;aCent[204][Y]=-4.999*SF;aCent[204][Z]=3.613*SF;aCent[204][W]=HYDR;
aCent[205][X]=8.121*SF;aCent[205][Y]=-6.159*SF;aCent[205][Z]=2.903*SF;aCent[205][W]=HYDR;
aCent[206][X]=5.393*SF;aCent[206][Y]=-7.382*SF;aCent[206][Z]=3.344*SF;aCent[206][W]=HYDR;
aCent[207][X]=6.880*SF;aCent[207][Y]=-8.302*SF;aCent[207][Z]=3.093*SF;aCent[207][W]=HYDR;
aCent[208][X]=7.707*SF;aCent[208][Y]=-7.394*SF;aCent[208][Z]=4.970*SF;aCent[208][W]=HYDR;
aCent[209][X]=6.637*SF;aCent[209][Y]=-5.290*SF;aCent[209][Z]=0.434*SF;aCent[209][W]=NITR;
aCent[210][X]=6.389*SF;aCent[210][Y]=-5.315*SF;aCent[210][Z]=-1.015*SF;aCent[210][W]=CARB;
aCent[211][X]=7.332*SF;aCent[211][Y]=-4.405*SF;aCent[211][Z]=-1.823*SF;aCent[211][W]=CARB;
aCent[212][X]=7.082*SF;aCent[212][Y]=-4.123*SF;aCent[212][Z]=-2.993*SF;aCent[212][W]=OXYG;
aCent[213][X]=4.914*SF;aCent[213][Y]=-4.993*SF;aCent[213][Z]=-1.265*SF;aCent[213][W]=CARB;
aCent[214][X]=4.431*SF;aCent[214][Y]=-5.743*SF;aCent[214][Z]=-2.358*SF;aCent[214][W]=OXYG;
aCent[215][X]=6.509*SF;aCent[215][Y]=-4.415*SF;aCent[215][Z]=0.930*SF;aCent[215][W]=HYDR;
aCent[216][X]=6.562*SF;aCent[216][Y]=-6.329*SF;aCent[216][Z]=-1.378*SF;aCent[216][W]=HYDR;
aCent[217][X]=4.344*SF;aCent[217][Y]=-5.236*SF;aCent[217][Z]=-0.389*SF;aCent[217][W]=HYDR;
aCent[218][X]=4.778*SF;aCent[218][Y]=-3.934*SF;aCent[218][Z]=-1.457*SF;aCent[218][W]=HYDR;
aCent[219][X]=3.714*SF;aCent[219][Y]=-6.324*SF;aCent[219][Z]=-1.987*SF;aCent[219][W]=HYDR;
aCent[220][X]=8.419*SF;aCent[220][Y]=-3.920*SF;aCent[220][Z]=-1.202*SF;aCent[220][W]=NITR;
aCent[221][X]=9.451*SF;aCent[221][Y]=-3.116*SF;aCent[221][Z]=-1.870*SF;aCent[221][W]=CARB;
aCent[222][X]=8.984*SF;aCent[222][Y]=-1.725*SF;aCent[222][Z]=-2.316*SF;aCent[222][W]=CARB;
aCent[223][X]=9.539*SF;aCent[223][Y]=-1.177*SF;aCent[223][Z]=-3.267*SF;aCent[223][W]=OXYG;
aCent[224][X]=8.573*SF;aCent[224][Y]=-4.210*SF;aCent[224][Z]=-0.246*SF;aCent[224][W]=HYDR;
aCent[225][X]=10.297*SF;aCent[225][Y]=-2.987*SF;aCent[225][Z]=-1.194*SF;aCent[225][W]=HYDR;
aCent[226][X]=9.805*SF;aCent[226][Y]=-3.652*SF;aCent[226][Z]=-2.752*SF;aCent[226][W]=HYDR;
aCent[227][X]=7.956*SF;aCent[227][Y]=-1.164*SF;aCent[227][Z]=-1.660*SF;aCent[227][W]=NITR;
aCent[228][X]=7.289*SF;aCent[228][Y]=0.084*SF;aCent[228][Z]=-2.054*SF;aCent[228][W]=CARB;
aCent[229][X]=6.855*SF;aCent[229][Y]=0.916*SF;aCent[229][Z]=-0.829*SF;aCent[229][W]=CARB;
aCent[230][X]=6.222*SF;aCent[230][Y]=0.366*SF;aCent[230][Z]=0.076*SF;aCent[230][W]=OXYG;
aCent[231][X]=6.110*SF;aCent[231][Y]=-0.243*SF;aCent[231][Z]=-2.994*SF;aCent[231][W]=CARB;
aCent[232][X]=5.046*SF;aCent[232][Y]=-1.171*SF;aCent[232][Z]=-2.378*SF;aCent[232][W]=CARB;
aCent[233][X]=3.923*SF;aCent[233][Y]=-1.592*SF;aCent[233][Z]=-3.338*SF;aCent[233][W]=CARB;
aCent[234][X]=4.251*SF;aCent[234][Y]=-2.811*SF;aCent[234][Z]=-4.100*SF;aCent[234][W]=NITR;
aCent[235][X]=4.859*SF;aCent[235][Y]=-2.914*SF;aCent[235][Z]=-5.274*SF;aCent[235][W]=CARB;
aCent[236][X]=5.289*SF;aCent[236][Y]=-1.864*SF;aCent[236][Z]=-5.937*SF;aCent[236][W]=NITR;
aCent[237][X]=5.035*SF;aCent[237][Y]=-4.095*SF;aCent[237][Z]=-5.809*SF;aCent[237][W]=NITR;
aCent[238][X]=7.579*SF;aCent[238][Y]=-1.676*SF;aCent[238][Z]=-0.874*SF;aCent[238][W]=HYDR;
aCent[239][X]=8.009*SF;aCent[239][Y]=0.663*SF;aCent[239][Z]=-2.630*SF;aCent[239][W]=HYDR;
aCent[240][X]=5.634*SF;aCent[240][Y]=0.678*SF;aCent[240][Z]=-3.269*SF;aCent[240][W]=HYDR;
aCent[241][X]=6.524*SF;aCent[241][Y]=-0.720*SF;aCent[241][Z]=-3.880*SF;aCent[241][W]=HYDR;
aCent[242][X]=5.538*SF;aCent[242][Y]=-2.059*SF;aCent[242][Z]=-2.031*SF;aCent[242][W]=HYDR;
aCent[243][X]=4.579*SF;aCent[243][Y]=-0.652*SF;aCent[243][Z]=-1.549*SF;aCent[243][W]=HYDR;
aCent[244][X]=3.033*SF;aCent[244][Y]=-1.774*SF;aCent[244][Z]=-2.766*SF;aCent[244][W]=HYDR;
aCent[245][X]=3.669*SF;aCent[245][Y]=-0.765*SF;aCent[245][Z]=-4.003*SF;aCent[245][W]=HYDR;
aCent[246][X]=3.963*SF;aCent[246][Y]=-3.694*SF;aCent[246][Z]=-3.698*SF;aCent[246][W]=HYDR;
aCent[247][X]=5.150*SF;aCent[247][Y]=-0.962*SF;aCent[247][Z]=-5.521*SF;aCent[247][W]=HYDR;
aCent[248][X]=5.761*SF;aCent[248][Y]=-1.962*SF;aCent[248][Z]=-6.815*SF;aCent[248][W]=HYDR;
aCent[249][X]=4.649*SF;aCent[249][Y]=-4.894*SF;aCent[249][Z]=-5.327*SF;aCent[249][W]=HYDR;
aCent[250][X]=5.508*SF;aCent[250][Y]=-4.205*SF;aCent[250][Z]=-6.684*SF;aCent[250][W]=HYDR;
aCent[251][X]=7.156*SF;aCent[251][Y]=2.230*SF;aCent[251][Z]=-0.780*SF;aCent[251][W]=NITR;
aCent[252][X]=6.782*SF;aCent[252][Y]=3.088*SF;aCent[252][Z]=0.345*SF;aCent[252][W]=CARB;
aCent[253][X]=5.261*SF;aCent[253][Y]=3.331*SF;aCent[253][Z]=0.395*SF;aCent[253][W]=CARB;
aCent[254][X]=4.586*SF;aCent[254][Y]=3.165*SF;aCent[254][Z]=-0.624*SF;aCent[254][W]=OXYG;
aCent[255][X]=7.554*SF;aCent[255][Y]=4.394*SF;aCent[255][Z]=0.119*SF;aCent[255][W]=CARB;
aCent[256][X]=7.677*SF;aCent[256][Y]=4.474*SF;aCent[256][Z]=-1.401*SF;aCent[256][W]=CARB;
aCent[257][X]=7.820*SF;aCent[257][Y]=3.010*SF;aCent[257][Z]=-1.816*SF;aCent[257][W]=CARB;
aCent[258][X]=7.107*SF;aCent[258][Y]=2.628*SF;aCent[258][Z]=1.279*SF;aCent[258][W]=HYDR;
aCent[259][X]=7.009*SF;aCent[259][Y]=5.234*SF;aCent[259][Z]=0.505*SF;aCent[259][W]=HYDR;
aCent[260][X]=8.548*SF;aCent[260][Y]=4.308*SF;aCent[260][Z]=0.561*SF;aCent[260][W]=HYDR;
aCent[261][X]=6.800*SF;aCent[261][Y]=4.914*SF;aCent[261][Z]=-1.836*SF;aCent[261][W]=HYDR;
aCent[262][X]=8.540*SF;aCent[262][Y]=5.066*SF;aCent[262][Z]=-1.707*SF;aCent[262][W]=HYDR;
aCent[263][X]=7.349*SF;aCent[263][Y]=2.844*SF;aCent[263][Z]=-2.766*SF;aCent[263][W]=HYDR;
aCent[264][X]=8.876*SF;aCent[264][Y]=2.739*SF;aCent[264][Z]=-1.855*SF;aCent[264][W]=HYDR;
aCent[265][X]=4.710*SF;aCent[265][Y]=3.739*SF;aCent[265][Z]=1.555*SF;aCent[265][W]=NITR;
aCent[266][X]=3.287*SF;aCent[266][Y]=4.031*SF;aCent[266][Z]=1.686*SF;aCent[266][W]=CARB;
aCent[267][X]=2.901*SF;aCent[267][Y]=5.305*SF;aCent[267][Z]=0.913*SF;aCent[267][W]=CARB;
aCent[268][X]=3.684*SF;aCent[268][Y]=6.256*SF;aCent[268][Z]=0.871*SF;aCent[268][W]=OXYG;
aCent[269][X]=3.035*SF;aCent[269][Y]=4.190*SF;aCent[269][Z]=3.187*SF;aCent[269][W]=CARB;
aCent[270][X]=4.385*SF;aCent[270][Y]=4.655*SF;aCent[270][Z]=3.729*SF;aCent[270][W]=CARB;
aCent[271][X]=5.393*SF;aCent[271][Y]=3.949*SF;aCent[271][Z]=2.823*SF;aCent[271][W]=CARB;
aCent[272][X]=2.719*SF;aCent[272][Y]=3.181*SF;aCent[272][Z]=1.316*SF;aCent[272][W]=HYDR;
aCent[273][X]=2.274*SF;aCent[273][Y]=4.924*SF;aCent[273][Z]=3.372*SF;aCent[273][W]=HYDR;
aCent[274][X]=2.781*SF;aCent[274][Y]=3.223*SF;aCent[274][Z]=3.618*SF;aCent[274][W]=HYDR;
aCent[275][X]=4.482*SF;aCent[275][Y]=5.721*SF;aCent[275][Z]=3.654*SF;aCent[275][W]=HYDR;
aCent[276][X]=4.518*SF;aCent[276][Y]=4.377*SF;aCent[276][Z]=4.775*SF;aCent[276][W]=HYDR;
aCent[277][X]=6.262*SF;aCent[277][Y]=4.562*SF;aCent[277][Z]=2.682*SF;aCent[277][W]=HYDR;
aCent[278][X]=5.662*SF;aCent[278][Y]=2.983*SF;aCent[278][Z]=3.253*SF;aCent[278][W]=HYDR;
aCent[279][X]=1.688*SF;aCent[279][Y]=5.360*SF;aCent[279][Z]=0.336*SF;aCent[279][W]=NITR;
aCent[280][X]=1.185*SF;aCent[280][Y]=6.543*SF;aCent[280][Z]=-0.353*SF;aCent[280][W]=CARB;
aCent[281][X]=0.715*SF;aCent[281][Y]=7.607*SF;aCent[281][Z]=0.655*SF;aCent[281][W]=CARB;
aCent[282][X]=1.513*SF;aCent[282][Y]=7.324*SF;aCent[282][Z]=-0.353*SF;aCent[282][W]=OXYG;
aCent[283][X]=0.048 *SF;aCent[283][Y]=6.014*SF;aCent[283][Z]=-1.229*SF;aCent[283][W]=CARB;
aCent[284][X]=-0.519 *SF;aCent[284][Y]=4.852*SF;aCent[284][Z]=-0.412*SF;aCent[284][W]=CARB;
aCent[285][X]=0.716*SF;aCent[285][Y]=4.275*SF;aCent[285][Z]=0.272*SF;aCent[285][W]=CARB;
aCent[286][X]=1.961*SF;aCent[286][Y]=6.966*SF;aCent[286][Z]=-0.991*SF;aCent[286][W]=HYDR;
aCent[287][X]= -0.697*SF;aCent[287][Y]=6.770*SF;aCent[287][Z]=-1.389*SF;aCent[287][W]=HYDR;
aCent[288][X]= 0.463*SF;aCent[288][Y]=5.630*SF;aCent[288][Z]=-2.162*SF;aCent[288][W]=HYDR;
aCent[289][X]=-1.232 *SF;aCent[289][Y]= 5.201*SF;aCent[289][Z]= 0.310*SF;aCent[289][W]=HYDR;
aCent[290][X]=-1.019*SF;aCent[290][Y]=4.114*SF;aCent[290][Z]=-1.041*SF;aCent[290][W]=HYDR;
aCent[291][X]=0.470*SF;aCent[291][Y]=3.937*SF;aCent[291][Z]=1.260*SF;aCent[291][W]=HYDR;
aCent[292][X]=1.121*SF;aCent[292][Y]=3.461*SF;aCent[292][Z]=-0.329 *SF;aCent[292][W]=HYDR;
aCent[293][X]=1.271*SF;aCent[293][Y]=8.822*SF;aCent[293][Z]=0.549*SF;aCent[293][W]=NITR;
aCent[294][X]=0.852*SF;aCent[294][Y]=10.027 *SF;aCent[294][Z]=1.285*SF;aCent[294][W]=CARB;
aCent[295][X]=-0.406*SF;aCent[295][Y]=10.657*SF;aCent[295][Z]=0.683*SF;aCent[295][W]=CARB;
aCent[296][X]=-0.387*SF;aCent[296][Y]=10.916*SF;aCent[296][Z]=-0.540*SF;aCent[296][W]=OXYG;
aCent[297][X]=1.972*SF;aCent[297][Y]=11.071*SF;aCent[297][Z]=1.284*SF;aCent[297][W]=CARB;
aCent[298][X]=3.120*SF;aCent[298][Y]=10.541*SF;aCent[298][Z]=1.911*SF;aCent[298][W]=OXYG;
aCent[299][X]=-1.341*SF;aCent[299][Y]=10.903*SF;aCent[299][Z]=1.473*SF;aCent[299][W]=OXYG;
aCent[300][X]=1.969*SF;aCent[300][Y]=8.961*SF;aCent[300][Z]=-0.165*SF;aCent[300][W]=HYDR;
aCent[301][X]=0.601*SF;aCent[301][Y]=9.760*SF;aCent[301][Z]=2.310*SF;aCent[301][W]=HYDR;
aCent[302][X]=2.210*SF;aCent[302][Y]=11.338*SF;aCent[302][Z]=0.272*SF;aCent[302][W]=HYDR;
aCent[303][X]=1.636*SF;aCent[303][Y]=11.959*SF;aCent[303][Z]=1.824*SF;aCent[303][W]=HYDR;
aCent[304][X]=2.831*SF;aCent[304][Y]=10.040*SF;aCent[304][Z]=2.676*SF;aCent[304][W]=HYDR;
/**/
SF=1;
	FOR(ii,1,305){
		atomInfoNum[ii][ATOTYP]=(int)aCent[ii][W];
	}
alphaCIndices[1]=1;		/*ASN*/  strcpy(atomInfoStr[1][RESTYP],"ASN"); residueType[1]=ASN;
alphaCIndices[2]=17;    /*LEU*/	strcpy(atomInfoStr[17][RESTYP],"LEU"); residueType[2]=LEU;
alphaCIndices[3]=36;    /*TYR*/strcpy(atomInfoStr[36][RESTYP],"TYR"); residueType[3]=TYR;
alphaCIndices[4]=57;    /*ILE*/strcpy(atomInfoStr[57][RESTYP],"ILE"); residueType[4]=ILE;
alphaCIndices[5]=76;    /*GLN*/strcpy(atomInfoStr[76][RESTYP],"GLN"); residueType[5]=GLN;
alphaCIndices[6]=93;    /*TRP*/strcpy(atomInfoStr[93][RESTYP],"TRP"); residueType[6]=TRP;
alphaCIndices[7]=117;   /*LEU*/strcpy(atomInfoStr[117][RESTYP],"LEU"); residueType[7]=LEU;
alphaCIndices[8]=136;   /* LYS*/strcpy(atomInfoStr[136][RESTYP],"LYS"); residueType[8]=LYS;
alphaCIndices[9]=158;   /* ASP*/strcpy(atomInfoStr[158][RESTYP],"ASP"); residueType[9]=ASP;
alphaCIndices[10]=170;   /* GLY*/strcpy(atomInfoStr[170][RESTYP],"GLY"); residueType[10]=GLY;
alphaCIndices[11]=177;  /*  GLY*/strcpy(atomInfoStr[177][RESTYP],"GLY"); residueType[11]=GLY;
alphaCIndices[12]=184;  /*  PRO*/strcpy(atomInfoStr[184][RESTYP],"PRO"); residueType[12]=PRO;
alphaCIndices[13]=198;  /*  SER*/strcpy(atomInfoStr[198][RESTYP],"SER"); residueType[13]=SER;
alphaCIndices[14]=209;  /*  SER*/strcpy(atomInfoStr[209][RESTYP],"SER"); residueType[14]=SER;
alphaCIndices[15]=220;  /*  GLY*/strcpy(atomInfoStr[220][RESTYP],"GLY"); residueType[15]=GLY;
alphaCIndices[16]=227;  /*  ARG*/strcpy(atomInfoStr[227][RESTYP],"ARG"); residueType[16]=ARG;
alphaCIndices[17]=251;  /*  PRO*/strcpy(atomInfoStr[251][RESTYP],"PRO"); residueType[17]=PRO;
alphaCIndices[18]=265;  /*  PRO*/strcpy(atomInfoStr[265][RESTYP],"PRO"); residueType[18]=PRO;
alphaCIndices[19]=279;  /*  PRO*/strcpy(atomInfoStr[279][RESTYP],"PRO"); residueType[19]=PRO;
alphaCIndices[20]=293;  /*  SER*/strcpy(atomInfoStr[293][RESTYP],"SER"); residueType[20]=SER;
/**/FOR(ii,0,20){
		alphaCIndices[ii+1]+=1;/*fix previous misnaming*/
		strcpy(atomInfoStr[alphaCIndices[ii+1]][RESTYP], atomInfoStr[alphaCIndices[ii+1]-1][RESTYP]);
		bbIndices[1][3*ii+1]=alphaCIndices[ii+1]-1;
		bbIndices[1][3*ii+2]=alphaCIndices[ii+1];
		bbIndices[1][3*ii+3]=alphaCIndices[ii+1]+1;


	}/**/
/**/setAtomInfoNumRESTYN();/**/
spineKnots=SPINENUM;
backBone=1;
flyForward=0;
/**/targetAC=6;/*TRP*/
lastAtom=304;
lastResidue=20;
lastChain=1;
loxLength=40;
sceneType=/**LOXODROMIC**/MOLECULE;
grids();
initializeFigures();
plane();
constructLox();

/**/ constructACircle();/**/
strcpy(showMolName , "Gila Monster Venom");
constructPrototypeAtomSphere();
/**constructMolecule();**/
setAtomRadii();
reConstructMarker();
lastBB[1]=20*3;
noConnections=0;
computeBackBone();
/**computeFrame();**/
reSet();
test100();
/*end akgeowall*/
}
void switchAffToAff2(void){
	int ii;
	FOR(ii,0,12){
		aff2[ii]=aff[ii];
	}aff2[12]=newCent[X];aff2[13]=newCent[Y];aff2[14]=newCent[Z];
}
void switchAff2ToAff(void){
	int ii;
	FOR(ii,0,12){
		aff[ii]=aff2[ii];
	}aff[12]=newCent[X];aff[13]=newCent[Y];aff[14]=newCent[Z];
}
/***der_Gamma is same for both Bishop and Frenet frames***/
float der_Gamma(float TTT){
        return  TTT;
}
/**********************************************************************/

/**********************************************************************/
float der_TttF(int ttt,float TTT,float NNN,float BBB){
        return NNN*development1[ttt][1];
         /* Normal*kappa(ttt) */
}
/**********************************************************************/
float der_NnnF(int ttt,float TTT,float NNN,float BBB){
        return  -development1[ttt][1]*TTT + development2[ttt][1]*BBB;
        /*return  -kappa(ttt)*TTT + tau(ttt)*BBB;*/

}
/**********************************************************************/
float der_BbbF(int ttt,float TTT,float NNN,float BBB){
        return -development2[ttt][1]*NNN;
        /*return -tau(ttt)*NNN;*/
}
/**********************************************************************/

/**********************************************************************/
float der_TttB(int ttt,float TTT,float NNN,float BBB){
    return development1[ttt][1]*NNN + development2[ttt][1]*BBB;
    /*return  kappa(ttt)*NNN + tau(ttt)*BBB;*/
}
/**********************************************************************/
float der_NnnB(int ttt,float TTT,float NNN,float BBB){
   return  -1*development1[ttt][1]*TTT;
   /*return  -kappa(ttt)*TTT;*/
}
/**********************************************************************/
float der_BbbB(int ttt,float TTT,float NNN,float BBB){
   return  -1*development2[ttt][1]*TTT;
   /*return  -tau(ttt)*TTT;*/
}
/***** end of McCreary's curve/frame computation functions *************/
void calculatefrenetcurve(void){ /* calculates the curve with Runge-Kutta */
   for(jj=0; jj<3; jj++){ /* for X, Y, and Z*/
      /**derivative functions for T, N, & B use precomputed values of k_1 and k_2**/
      /**to construct the Bishop or Frenet framed curve**/
      /**Note: k_1=kappa; k_2=tau if whichFrame==Frenet**/
      int ii;

      /**(G)amma curve**/
      float  G1, dG, g1, g2, g3, g4;
      /** initial value, step, variables for transition functions:**/
      /** 4 transition variables for this Runge-Kutta **/
      /**Bishop frame vectors (T,N,B). (n.b., need all three)**/
      float  T1, dT, tmpT, t1, t2, t3, t4,
             N1, dN, tmpN, n1, n2, n3, n4,
             B1, dB, tmpB, b1, b2, b3, b4;
      /**time variables**/
      float  h,hhh;
      /**parameter variables**/
      int  tt1,tmp1;

      /**Following is (potentially variable) number of points placed between each point computed on Gamma**/
      /**As ppi gets larger, this is equivalent to decreasing the step size in the R-K process.**/
      /**Note that the same number of points on the curve will be computed no matter what ppi's value **/
      int ppi=1;

      /**** have initframe initialized in deFault or somewhere else**/
      /*hhh=portion of curve traversed each step*/
      hhh=(float)1.0/(numbpts-1); /*gammaLength;*/
      h=fabs(hhh);/**only positive steps here. Safety check**/

      h*= paramlength;  /* scale the stepsize */

      /** initial frame values specified in determineInitialConditions() **/
      tt1=0;
      G1=initpos[jj];
      T1=initframe[jj+0];
      N1=initframe[jj+3];
      B1=initframe[jj+6];

      FOR(ii,0,numbpts*ppi){
         /*record the solution for this step**/
         if(ii%ppi==0){
            /* the frame points are located in space centered at point on curve, not centered at origin*/
            /* gammaCurve[whichFrame][ii][whichCoord]=G1; */
            /* frame[whichFrame][TT][ii/ppi][whichCoord]=T1+gammaCurve[whichFrame][ii][whichCoord]; */
            /* frame[whichFrame][NN][ii/ppi][whichCoord]=N1+gammaCurve[whichFrame][ii][whichCoord]; */
            /* frame[whichFrame][BB][ii/ppi][whichCoord]=B1+gammaCurve[whichFrame][ii][whichCoord]; */

            /*   frenetframe[0 3 6]  [T0 N0 B0] */
            /*              [1 4 7]==[T1 N1 B1] */
            /*              [2 5 8]  [T2 N2 B2] */

            frenetcurve[ii/ppi][jj]        = G1;  /* Scale things to be huge! */
            curvefrenetframe[ii/ppi][jj+0] = T1;
            curvefrenetframe[ii/ppi][jj+3] = N1;
            curvefrenetframe[ii/ppi][jj+6] = B1;

         }
         /*step 1*/
         tt1 = 2*ii;
         t1=h*der_TttF(tt1,T1,N1,B1);
         n1=h*der_NnnF(tt1,T1,N1,B1);
         b1=h*der_BbbF(tt1,T1,N1,B1);
         g1=h*der_Gamma(T1);

         /*step 2*/
         tmp1=tt1+1;
           tmpT=T1+.5*t1;
           tmpN=N1+.5*n1;
           tmpB=B1+.5*b1;

         t2=h*der_TttF(tmp1,tmpT,tmpN,tmpB);
         n2=h*der_NnnF(tmp1,tmpT,tmpN,tmpB);
         b2=h*der_BbbF(tmp1,tmpT,tmpN,tmpB);
         g2=h*der_Gamma(tmpT);

         /*step 3*/
           tmpT=T1+.5*t2;
           tmpN=N1+.5*n2;
           tmpB=B1+.5*b2;
         t3=h*der_TttF(tmp1,tmpT,tmpN,tmpB);
         n3=h*der_NnnF(tmp1,tmpT,tmpN,tmpB);
         b3=h*der_BbbF(tmp1,tmpT,tmpN,tmpB);
         g3=h*der_Gamma(tmpT);

         /*step 4*/
         tmp1=tt1+2;
           tmpT=T1+t3;
           tmpN=N1+n3;
           tmpB=B1+b3;
         t4=h*der_TttF(tmp1,tmpT,tmpN,tmpB);
         n4=h*der_NnnF(tmp1,tmpT,tmpN,tmpB);
         b4=h*der_BbbF(tmp1,tmpT,tmpN,tmpB);
         g4=h*der_Gamma(tmpT);

         /*changes in all parameters*/
         dG=(g1+g2+g2+g3+g3+g4)/6.;
         dT=(t1+t2+t2+t3+t3+t4)/6.;
         dN=(n1+n2+n2+n3+n3+n4)/6.;
         dB=(b1+b2+b2+b3+b3+b4)/6.;

         /*update parameter values*/
         G1+=dG;
         T1+=dT;
         N1+=dN;
         B1+=dB;
     }/*end of loop*/
   }
}
/**********************************************************************/
/* calculatebishopcurve - calculates the curve if the developments are*/
/*    the Bishop parameters.                                          */
/**********************************************************************/
void calculatebishopcurve(void){ /* calculates the curve with Runge-Kutta */
      /**derivative functions for T, N, & B use precomputed values of k_1 and k_2**/
      /**to construct the Bishop or Frenet framed curve**/
      /**Note: k_1=kappa; k_2=tau if whichFrame==Frenet**/
      int ii;

      /**(G)amma curve**/
      float  G1, dG, g1, g2, g3, g4;
      /** initial value, step, variables for transition functions:**/
      /** 4 transition variables for this Runge-Kutta **/
      /**Bishop frame vectors (T,N,B). (n.b., need all three)**/
      float  T1, dT, tmpT, t1, t2, t3, t4,
             N1, dN, tmpN, n1, n2, n3, n4,
             B1, dB, tmpB, b1, b2, b3, b4;
      /**time variables**/
      float  h,hhh;
      /**parameter variables**/
      int  tt1,tmp1;

      /**Following is (potentially variable) number of points placed between each point computed on Gamma**/
      /**As ppi gets larger, this is equivalent to decreasing the step size in the R-K process.**/
      /**Note that the same number of points on the curve will be computed no matter what ppi's value **/
      int ppi=1;

      /**** have initframe initialized in deFault or somewhere else**/
      /*hhh=portion of curve traversed each step*/
      hhh=(float)1.0/(numbpts-1); /*gammaLength;*/
      h=fabs(hhh);/**only positive steps here. Safety check**/

      h*= paramlength;  /* scale the stepsize */


   for(jj=0; jj<3; jj++){ /* for X, Y, and Z*/

      /** initial frame values specified in determineInitialConditions() **/
      tt1=0;
      G1=initpos[jj];
      T1=initframe[jj+0];
      N1=initframe[jj+3];
      B1=initframe[jj+6];

      FOR(ii,0,numbpts*ppi){
         /*record the solution for this step**/
         if(ii%ppi==0){
            /* the frame points are located in space centered at point on curve, not centered at origin*/
            /* gammaCurve[whichFrame][ii][whichCoord]=G1; */
            /* frame[whichFrame][TT][ii/ppi][whichCoord]=T1+gammaCurve[whichFrame][ii][whichCoord]; */
            /* frame[whichFrame][NN][ii/ppi][whichCoord]=N1+gammaCurve[whichFrame][ii][whichCoord]; */
            /* frame[whichFrame][BB][ii/ppi][whichCoord]=B1+gammaCurve[whichFrame][ii][whichCoord]; */

            /*   showbishopframe[0 3 6]  [T0 U0 V0] */
            /*              [1 4 7]==[T1 U1 V1] */
            /*              [2 5 8]  [T2 U2 V2] */

            bishopcurve[ii/ppi][jj]        = G1;
            curvebishopframe[ii/ppi][jj+0] = T1;
            curvebishopframe[ii/ppi][jj+3] = N1;
            curvebishopframe[ii/ppi][jj+6] = B1;

         }
         /*step 1*/
         tt1 = 2*ii;
         t1=h*der_TttB(tt1,T1,N1,B1);
         n1=h*der_NnnB(tt1,T1,N1,B1);
         b1=h*der_BbbB(tt1,T1,N1,B1);
         g1=h*der_Gamma(T1);

         /*step 2*/
         tmp1=tt1+1;
           tmpT=T1+.5*t1;
           tmpN=N1+.5*n1;
           tmpB=B1+.5*b1;

         t2=h*der_TttB(tmp1,tmpT,tmpN,tmpB);
         n2=h*der_NnnB(tmp1,tmpT,tmpN,tmpB);
         b2=h*der_BbbB(tmp1,tmpT,tmpN,tmpB);
         g2=h*der_Gamma(tmpT);

         /*step 3*/
           tmpT=T1+.5*t2;
           tmpN=N1+.5*n2;
           tmpB=B1+.5*b2;
         t3=h*der_TttB(tmp1,tmpT,tmpN,tmpB);
         n3=h*der_NnnB(tmp1,tmpT,tmpN,tmpB);
         b3=h*der_BbbB(tmp1,tmpT,tmpN,tmpB);
         g3=h*der_Gamma(tmpT);

         /*step 4*/
         tmp1=tt1+2;
           tmpT=T1+t3;
           tmpN=N1+n3;
           tmpB=B1+b3;
         t4=h*der_TttB(tmp1,tmpT,tmpN,tmpB);
         n4=h*der_NnnB(tmp1,tmpT,tmpN,tmpB);
         b4=h*der_BbbB(tmp1,tmpT,tmpN,tmpB);
         g4=h*der_Gamma(tmpT);

         /*changes in all parameters*/
         dG=(g1+g2+g2+g3+g3+g4)/6.;
         dT=(t1+t2+t2+t3+t3+t4)/6.;
         dN=(n1+n2+n2+n3+n3+n4)/6.;
         dB=(b1+b2+b2+b3+b3+b4)/6.;

         /*update parameter values*/
         G1+=dG;
         T1+=dT;
         N1+=dN;
         B1+=dB;
     }/*end of loop*/
     /*******constructCurve();/*and frame from the frame and path coordinates*/
   }
}



/**********************************************************************/
/* Code for stars - do not meddle in the afairs of wizards...         */
/**********************************************************************/

#ifdef IRIX
#include <alloca.h>
#endif

#ifndef OSXMac

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <GL/gl.h>

#define NewA(type, count)  ((type *)alloca((count)*sizeof(type)))
#define NewN(type, count)  ((type *)malloc((count)*sizeof(type)))

typedef float Matrix[4*4];
typedef struct { float x[3]; } Point;

int zoom2 = 0;		/* Flag for "enlarge stars twofold" */

 extern int    caveyes;	    /* Hope this is defined somewhere else*/
 extern Matrix s_starmat;  /* Hope this is defined elsewhere too */

void vtfmvector( Point *dst, register Point *src, register Matrix T )
{
  int i;
  for(i = 0; i < 3; i++)
    dst->x[i] = src->x[0]*T[i] + src->x[1]*T[i+4] + src->x[2]*T[i+8];
}

struct starsize {
  float size;
  int base;
  int count;
  int color;
} *starsizes;
int nstarsizes = 0;

float scolors[][3] = {
	.5, 1, 1,	/* W */
	.7, .8, 1,	/* O */
	.9, .9, 1,	/* B */
	1, 1, 1,	/* A */
	1, 1, .9,	/* F */
	1, 1, .8,	/* G */
	1, .8, .7,	/* K */
	1, .7, .5,	/* M */
};

float *starpos;

#define SRADIX  91
#define SZERO   '$'
#define SSIZE   '!'
#define SSIZESCL 8
#define SRASCALE (8192 / (2*M_PI))
#define SDECSCALE (8192 / M_PI)
#define SDECBIAS 4096
#define SCOLORBASE 0
#define SSIZEBASE  8

#define SRADIX  91
#define SZERO   '$'
#define SSIZE   '!'
#define SSIZESCL 8
#define SRASCALE (8192 / (2*M_PI))
#define SDECSCALE (8192 / M_PI)
#define SDECBIAS 4096
#define SCOLORBASE 0
#define SSIZEBASE  8


/* 1128 real stars, to about visual magnitude 4.7 */
/* from  star2illi -m 4.7 -d 1.2 -c 3.5 */
char stardata[] =
  /* Show Pleiades at upper right, Orion at center of field */
"#@ -.98525 .062202 .1594 0  .16485 .095488 .98169 0  .045842 .99349 -.10433  0 0 0 1\n"
"!<'=BH_"
"!:(<$6`"
"!9'iodI!9&7`Lw!9)Z~2Y7mh&!9*YQZ["
"!8&*.4G!8(@eS\\!8+:7Te"
"!7'nLUM!7&Xh2oVCKK!7*5;Y=!7+arCl"
"!6'z/B6q[g`!6&S~33J'W$!6*A,_'Z~2Y"
"!5'@J`{@J`{FY.1!5&>3BTRa1Mes>N8BT58H_A9&PIw'9S!5+R{4M"
"!4'TGm$:OgQ<sY7Dm5_ezW?!4&i&?nWgiap]4^;xH'j|Cr$T_W9gL4Ra1M!4(0jiz>kCn"
"f-;R-Q}`!4*MQouCJ3:c).PGPLb+z\\h&fH&Xy>o[bv-!4+(Ebo!4%9?P'!4$B[9C"
"!3'PBX?S[8S^K^FFo3EV;lPMFm7!3&/ieQ8iPrW;6=[09AZg;y`'EdfJ=Q!3(pIe1$XnZ!3*"
"F<;>&Tm>g@ji+hf5c2?r%_;ruTU{!3+y59W!3%B9=%"
"!2'dJI2P]kssyp@N1[=bCKfkKB*:Pc]+3[J)Eo0!2&?iBE'VoFzXX\\G15QX99ER$H9QU7_"
"]BLA`DG.99?}\\4;LSC.Lec>D!2(8mH-!2)S3E@X4Z7LK8?!2*qwb$JiZy_'T96XaZ?B>N"
"[B^VhoB,n2VA!2+zR_)/HS)"
"!1'V)>_[cH~7:NOueHz]M.`TSd4Hd0Ok[Wz,-bR!1&L60m^N<JQ{3]b;Bv$oX[2^`{28]*"
"n+gY2se&_zCzaCD7?{U2RtHhei8*9.[Z!1(BPDr_e1>+K28ktFR6tfx!1)aSokav[ieak3"
"%[*FbZ`nTvVQ8SF]tgN9/Vkjy7`/vtPr=;]YH]\\vUzEK!1*fPS?i<D?eD5:wb2s=]7cgtAm"
"Q_Ed!1+`wO3;z\\=2tJ;!1%8|N)"
"!0'ZSd3]WtyEdi(dd]K[30S!0&1yhwv1>A;kA}>LE,Ws;jSv.{c;=}OQ1Q\\7;{73e^<s;I"
"dAqtjFCS][<bofP[tWt@k6aD!0(fg<ypBI\\Gnjr6/TQ!0)l)rrEOS~sY`0a3N`!0*Msg<"
"m5_$L\\HwqE9E9z?/G+b8Ge4Pcb5%ddcK7-Eoc\\Udft>P|gwp@,;DhqOW!0+,eOSJvejh\\>]"
"27+u"
"!/'505P=S2'yyI-Q}mUN2XfVyPe4lXz/1<rqo/vJYfO.7R]yW7D!/&e1DR7WHwJE-~KR23"
"eD4n8>On]b:`+2pry.VKWr<3`&=mFL3Sk`NXj_ad3(W;Am6S!/(dT;G=CWNlrRX+*_mJV\\f!/)(09_ACDX3v1jCvoFDyT8]'6{fa^s'*mx]9aay_]A*RI("
"\\Ge7!/*&N`MXwCa]onRkdC4gNL0w6n/JY2Cqqoy(>Kx>EC(NJaW!/+;L\\<\\SDETPRe/XdK"
")T;DnzZidbX7"
"!.'HBU|I|YH_EO?w2T.9cI[?GY==mb$E~hZP./^Xlq7_2Xf_LS;[IQ{!.&i6:&IA5fYf9}"
",U7:4/@.c>>$7rMYDT6Vee2`zHf5.]^_2B](&@k|^bB/gy7{Dg@J67Sn2)]+6HR9!.(?S\\$"
"f.I@hpuG9UEkGb<jPKQuqEX@S\\P>GfpVS\\P>^$_Wu=Ha:>Ixscd(!.)4lZ[bcdOOAA+NOI\\"
"oW/wA']7p1Jg0mUT)eXcHY1i+;76gF_^nfT7!.*k/FNC(U[y]r/1{L0%8L[]]>s*0iBZR`6"
"^YBzRC2m);LxA,>&J4JpTv-9cK;of]0_4:XorM3kzrFK{CR_P1hw:Ol2s+f}!.+0OF0!.%"
"8{U}"
"!-'g}UlaGZZfiREr$L<xOj2m/jtFycJazR$wsPA9e7OJK;zSE8g'Vd<R;PbC\\O*36T$Du9}!-&0{UtqQX|f9h&2&a2XI;{g~_HEM2b^y^3?5C]KWU`hFFO2/]6S(,yS-swXK:[]%8^F`R3_oBHNX5ia;h4"
"Ri7u!-(xCn98r1iv{]c2+f@>S[?HznU^LV=^LV=K;3^DV9bKq8u_lXq0%BS[7N4nVQS!-)"
"tLEmkWF0lCkd?s^wJiZy9|FX59Af|\\h;n=t2;u@>`}13Vd=@b.)VMl3S^PI\\!-**|Kq1@L>"
"AS<e6weWF&9Mg.mM4LYlc03PuE*@C\\/y>m-hLr3YoshG/qgLK1HZLqb/[R)O4lY$i$[v2%0["
"ieLu3t;rg>c^5HIr8d?==uDw$Z:0I)^%/'LWnlbWzx:GH;PL!-+ynM8.al|n9Z=O8sb!-%"
"94Oe"
"!,'(=5Gck`PsiS^5=OC@pL9l]H,%^;3DmZ,NiV=rZe[)d8O.ukH<mGG@s,dV;lO_M@5yK\\m4R@$"
":M;Z=\\@hA$BTQg6o`LFao),U9}dZDS?DEmewFN1qgWROli<cxD;;{dF{+eu9;KMt?@/'o~ht{E3u"
"+RFN.s8*Ew/mK'+~[0>.,YazE3C2VClR`0n?`jG=j7->lb:jpr`6y\\J7}}TM5qOGDx_HQUD^a;7y"
"k:XVkFN1tugnxfP{6|o:P0T=TK4JgrR<q-V`2ybw3Vht3hMZR31$j}f~/oioNYT'Ym>*Z,jzZt8>"
"[:?K\\n:E&{]1@R^Ncg6J*PjE-|Q46TWj=uJ~C3*UDk;bWrXw]E3[s1HHGw3GM:GrNnH3YSN%t2Zw"
"uPnI"
"!+'*)ee,r.a/IE69*U__@Z+-}.s.(=+15WPAB9oKf)dP32IR./&^W/qkr>)ku=Cr;CVxG;0.<i^"
":`Up>QI6H}IYQCUFzB6]iA8Onf<+uW]pwScu|jSo3yi88;M.DE;S^A`dqvDF|lg5<A[/?$PiQL0c"
"_]H^_s^MyH(@&*pP:i\\_DESrb>?GwXM/2C>veJDwy*CP22\\~@+`vYUh)Wr?kYAKtjF[>wNmUyLW-"
"8bS~AnEXb1f91*n~9?P'K;cEz&@g{<LO4Z\\3hntbjNN^jk1wpspSqs`FtJ0A{3M~%;0W.`d41Ui/"
"25EGW\\@P`K>[a}?`buz'nKaOsrdeuciaAZ8{bd*:p/Ji"
"!*'(?|1(Ch^-592.<J+4%7<59V-5?e_CLf[E`VyK7AOZ9vs*ZUY:Rg}R2)CS/ecUR^{X\\Ql`WgP"
"fDJX.@V*0Q;O1EF55GW=7BL^B[9BTD<wU68)]J9*qyY*vxJ+'yT{-<U:1WQ74_\\J5_\\P@|_NK|57"
"PX@)a{Hdt3HZzd;;|&T73Jj53}UN4ZY}6Gr67}JJH+yceOB+f)=cle:[u\\oXviqCw1a[xQ@p|`f_"
"2+];DeReOXP[jxcNorm?'`BBGg\\QS%Hw]lcd5TG4@&UPDwJ9PE1.Z6'3[g;7_&wx`OFXqJPN(LlZ"
"+]RH<QE?DV3/Q61CTB8TYk4n^^;fe*JYhic(BVOSFj4:X\\:7]FAyeKS+uVYcuY@Q/}Zt9V00;O_j"
">,n9BR.cC?*<H*8BZ2:Ah)13&@as&e4=13i&3?\\)63UN7YJU:D?EE'SxRc_1Xu<Jek|Ah*[Jh|2;"
"jAcnjMZ.l>d+mWPFxEhs%|1S&I_a,AUL;U?L>%HPJ$V%mjZ&mpYio$?E|B>-"
"!)'8jZAUGN:cHV-6@N><c6SB2R4M=;ue+FMl0dZoXwtt%6={GLJ{vFb=]66?KDTBXGG\\LR)vG5R"
"{r\\e$/N$+.9r8yUi:2[2>{`1Av8GNH`k[`C&ec^+gJ`.zjve{M@g'X\\f/1<r20J{:}XHMI[.Y{=;"
"ZV7kadGkr$NQxAQ$'(Tm3{M5<tGu[,Wt_3TcpRa.qgXVsubNw~k0$wGU*CSi4'3D;~S@Bu<pDKRb"
"k+/Gm*]D5U<+7TK+@GErBe=6C3>aEi<GKg49LF(gb&F=g]Q`l:ucmHTd.X@m7YRLG?^.PN:JZb_u"
"vv=9y)g2ym@W.AGf0T_U2Dqk6Yku6`Qs=QR8N?O5QS7aae99zcE0;;X/;gnT<<@eF+,dF@jnI}b^"
"M~EZmVj/x90S"
"!('4_Xo=;W]?(9]AV=XDx4\\J$Po\\,NpsaV%wGdt|zIg<WTb@bC[A<D(Gyk(Jn4}L~]GWeYgaSF~"
"s8DRw(@S|vQv~&07(P`).Q_^3)1x5-B/9s4xGY?']f3Co6r{qj7(qrmnsFKBtWF.xKfX%1cI-Are"
"9od[H82CYdJB^rZqe:cZnoCKp$^xzdUf{)ie'.eW.fbV:ROO@Z?TP<[/V42SV<0gh8:'q}c;zSRx"
"%}1S&li2+]tO1/nN7uad<cET@/W%YDjw\\U^Q]&G/]e>X^R:TeKNOfgC,}amj$9HC7'?<><wSGvc8"
"[fZWa@Dxs5hq|,WH2x2@F3d:GwsyX$@RZ1:HuGtb"
"!''97S*GiPJJg5PL92fa[X'h&Bkt^\\pwph>x)ii}BBz1X<s:,^oF<D+H~9ffGu-n[])o5C2&Y9|"
"1y>D8KRW8|NZ@}BmI%l)_\\DC`kC(e~9ijy\\Ckl<ilWqtr=^WxAQ$<@MQDSq5GYO]IoJOM)>LN(2q"
"UT3([CYQ\\.q~]q=^mWDX{W\\t$BN3%*dE41j3B$GnMa1ll]I'q1b^+-R[8hME@bC[DgMHFh>?OK5s"
"TK3K^HK}_Wf:f0L{gVO3k'S.'Un[/7[b/7[b/{do2Itb:h[+AX86MXTbN7\\W_E^(`94.ahLoxrjl"
"!&'6z[m??C*?xi]@OBoSR=%^p:a_MD0gUS6g}UFiO1ri\\;qo,^u(\\]@2\\Dd5/XL6LV,:01N=xW["
"@4ATSw4iWb7>W}qF`88Gc0K\\dsc`h_CPzTj((N[T.8^s3$E$=*U|?1Cp?H>^De[hK845Ky3JQ*TA"
"SNM%Sa8ZW@lEjenelxQ4x-Qd}t]Y*qKa1,ik:uIUFLpk^h?mgSYGrq@tuIj["
"!%'7(Xe=}FzB,P@H3SCa>)_i$LUm>srq8XCu,G=+5/4+kBE4aXF5NX|ET6_G7BZGNEpSs2SViMu"
"gpB8mf`,sJV+u0M,vfOyvf4\\y:GZzuEk"
;

static int slow_stars = 0;

void initstars() {
  char *cp;
  float ra, dec, cd;
  int maxstars = 5000;
  int maxsizes = 32;
  float *stars = NewA(float, maxstars*3);
  float *starp;
  int nstars = 0;
  struct starsize *sizes = NewA(struct starsize, maxsizes);
  struct starsize *sizep;
  int i, nsizes = 0;
  char line[512];
  FILE *inf = fopen("stars.illi", "r");
  float r = caveyes ? 90. : 1.;
  const GLubyte *rend = glGetString(GL_RENDERER);
  Point sp;
  static Matrix startfm = { 0,0,-1,0, -1,0,0,0, 0,1,0,0, 0,0,0,1 };

  if(rend==NULL || (rend[0]=='X'||rend[0]=='G'))
    slow_stars = 1;

  starp = &stars[0];
  sizep = &sizes[0];
  sizep->size = 1;
  sizep->base = 0;
  sizep->color = 3;
  nsizes = 0;
  cp = stardata;
  do {
    if(inf) {
	line[sizeof(line)-1] = '\0';
	if(fgets(line, sizeof(line)-1, inf) == NULL)
	    break;
	cp = line;
    }
    for(;;) {
      if(*cp >= '$') {
	if(nstars >= maxstars-1) {
	    maxstars *= 3;
	    starp = NewA(float, 3*maxstars);
	    memcpy(starp, stars, nstars*3*sizeof(float));
	    stars = starp;
	    starp = &stars[nstars*3];
	}
	ra = ((cp[0] - SZERO)*SRADIX + (cp[1] - SZERO)) / SRASCALE;
	dec = ((cp[2] - SZERO)*SRADIX + (cp[3] - SZERO) - SDECBIAS) / SDECSCALE;
	cd = cos(dec);
	sp.x[0] /*starp[2]*/ = -r * cd * cos(ra);
	sp.x[1] /*starp[0]*/ = -r * cd * sin(ra);
	sp.x[2] /*starp[1]*/ = r * sin(dec);
	vtfmvector( (Point *)starp, &sp, startfm );

	cp += 4;
	starp += 3;
	nstars++;
      } else if(*cp == '#') {
	/* Comments.
	 * Also, "#@" hack allows specifying star transformation.
	 */
	if(cp[1] == '@') {
	  cp += 2;
	  for(i = 0; i<16; i++)
	    startfm[i] = strtod(cp, &cp);
	  cp--;
	}
	while(*++cp != '\n' && *cp != '\0')
	    ;
      } else if(*cp == '!') {
	sizep->count = nstars - sizep->base;
	if(sizep->count > 0)
	    nsizes++, sizep++;
	if(nsizes >= maxsizes) {
	    maxsizes *= 3;
	    sizep = NewA(struct starsize, maxsizes);
	    memcpy(sizep, sizes, nsizes*sizeof(struct starsize));
	    sizes = sizep;
	    sizep = &sizes[nsizes];
	}
	sizep->size = (cp[1] - SZERO) / (float)SSIZESCL;
	sizep->color = (cp[2] - SZERO);
	sizep->base = nstars;
	cp += 3;
      } else if(*cp == '\0') {
	break;
      } else {
	cp++;
      }
    }
  } while(inf || *cp != '\0');
  sizep->count = nstars - sizep->base;
  nstarsizes = nsizes+1;
  starsizes = NewN(struct starsize, nstarsizes);
  memcpy(starsizes, sizes, nstarsizes*sizeof(struct starsize));

  starpos = NewN(float, 3*nstars);
  memcpy(starpos, stars, nstars*3*sizeof(float));
}

  void drawstars(float * lstarmat) {
  int s, i;
  float v, size;
  struct starsize *sp;
  float *starp;

  if(getenv("NOSTARS"))
    return;
  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  if(!slow_stars || zoom2)
    glEnable(GL_POINT_SMOOTH);
  glDisable(GL_DEPTH_TEST);
  glMultMatrixf(lstarmat);
  for(i = nstarsizes; --i >= 0; ) {
    sp = &starsizes[i];
    size = zoom2 ? 2*sp->size : sp->size;
    glPointSize( size );
    v = size < 1 ? size : 1 - (ceil(size) - size) * .25;
    glColor3f(	v*scolors[sp->color][0],
		v*scolors[sp->color][1],
		v*scolors[sp->color][2] );
    glBegin(GL_POINTS);
    for(s = sp->count, starp = &starpos[3*sp->base]; --s >= 0; starp += 3)
	glVertex3fv( starp );
    glEnd();
  }
  glPopMatrix();
  glDisable(GL_POINT_SMOOTH);
  glEnable(GL_DEPTH_TEST);
}
#endif  /**ending the if-not-define-OSXMac clause*******************/
/**********************************************************************/
/* arguments - Handles the command line arguments                     */
/**********************************************************************/
#ifdef IRIX
void arguments(int argc,char **argv){
  extern char *optarg;
  extern int optind;
  int chi;
  /* w: needs ONE number after -w, c<nothing> means NO number follows*/
  while ((chi = getopt(argc,argv,"cw:d:g:")) != -1)
  switch(chi)  {
                case 'c': caveyes=1; break;
                case 'w': win=atoi(optarg); break;
                case 'g': gap0=atof(optarg); break;
                }
  if (optind!=argc) fprintf(stderr,"%s: Incorrect usage\n",argv[0]);
}
#endif
//#ifdef OSXMac
#ifdef LINUX
void arguments(int argc,char **argv){
  extern char *optarg;
  extern int optind;
  int chi;
  /* w: needs ONE number after -w, c<nothing> means NO number follows*/
  while ((chi = getopt(argc,argv,"cw:d:g:")) != -1)
  switch(chi)  {
                case 'c': caveyes=1; break;
                case 'w': win=atoi(optarg); break;
                case 'g': gap0=atof(optarg); break;
                }
  if (optind!=argc) fprintf(stderr,"%s: Incorrect usage\n",argv[0]);
}
#endif
/**********************************************************************/

enum{console,cave,cube} eyes;

#ifdef WIN32
void arguments(int argc,char **argv){           /* Pat Hanrahan 1989 */

     //     eyes=cube; //default for syzygy, don't even think of compiling in caveOGL
     //   eyes=console;  //fix this by fixing arguments
     //whoa .. this should not have worked since WIN32 is NOT defined
      while(--argc){++argv; if(argv[0][0]=='-')switch(argv[0][1]){
      case 'w': win =atoi(argv[1]); argv++; argc--; break;
      case 'c': eyes=console;break; 
	  case 'g': gap0 =atof(argv[1]); argv++; argc--; break;
      case 'L': lux[0]=atof(argv[1]);
                lux[1]=atof(argv[2]);
                lux[2]=atof(argv[3]); argv +=3; argc -=3; break;
   }}}
#endif

void centerTheModel(void){
	int ii;
	float xx=0, yy=0, zz=0;
	FOR(ii,1,lastAtom+1){/*compute center of molecule*/
		xx+=aCent[ii][X];
		yy+=aCent[ii][Y];
		zz+=aCent[ii][Z];
	}
	xx/=lastAtom;
	yy/=lastAtom;
	zz/=lastAtom;
	FOR(ii,1,lastAtom+1){
		aCent[ii][X]-=xx;
		aCent[ii][Y]-=yy;
		aCent[ii][Z]-=zz;
	}
}
/**********************************************************/
/*************************************************************/
/**************************************************************/



/**********************************************************/
/*************************************************************/
//can't seem to get this to work. I tried putting a valid DATAPATH on the shellpath,
//but maybe this beast looks into some other path tree
//Damn ... it's a hardwired macro here !!!!
//can't make DATAPATH work. Checked with concrete and that source has NO DATAPATH
//prefixed ... and thus has to use local place for files (for now)
//took out the DATAPATH
/**************************************************************/
void openpdbfile(void){
	FILE *openptr;
	int ii=0, jj=1, kk=1, chaiNum=0, aCInx=1;/*ii=line#, jj=bbAtom#w/iChain**/
	char firstField[80],blank[1024],serialNum[80],atomName[80],altLoc[80];
	char resName[80],chainID[80],resNum[80],iCode[80],coord[80];
	char occup[80],tempf[80],segID[80],elementName[80],charge[80];
	char lastChainID[80];
	int virg=1;

	nnnn=0;/*qwert*/

	lastResidue=0;

	switch(openFileNum){
//gave up on an actual date-file path, just cd into data/ and execute binary
		case 0: openptr = fopen("pdb1a91.ent", "r");
			strcpy(showMolName , "subunit C");
			aff[14]= -50.0; noConnections=0;
			break;
		case 1: openptr = fopen("pdb1abv.ent", "r");
			strcpy(showMolName , "N-term delta");
			aff[14]= -50.0; noConnections=0;
			break;
		case 2: openptr = fopen("pdb1b9u.ent", "r");
			strcpy(showMolName , "subunit. B");
			aff[14]= -50.0; noConnections=0;
			break;
		case 3: openptr = fopen("pdb1qo1.ent", "r");
			strcpy(showMolName , "yeast rotor motor");
			aff[14]= -110.0; speed=100; targetAC=490; noConnections=1;
			break;
		case 4: openptr = fopen("pdb1fx0.ent", "r"); /**Chain**/
					strcpy(showMolName , "spinach F1");
					aff[14]= -100.0;
					break;
//we don't use these because they crash the VE at Beckman
		case 5: openptr = fopen("1I6X.pdb", "r"); /**Chain**/
					strcpy(showMolName , "CRP");
					aff[14]= -30.0;noConnections=0;
					break;
		case 6: openptr = fopen("cap.pdb", "r"); /**Chain**/
					strcpy(showMolName , "CRP from VMD");
					aff[14]= -40.0;noConnections=0;
					break;
		case 7: openptr = fopen("x.pdb", "r"); /**Chain**/
					strcpy(showMolName , "x.pdb");
					aff[14]= -40.0;
					break;

		case 10: openptr = fopen("pdb1a91.ent", "r");
			strcpy(showMolName , "subunit C");
			aff[14]= -50.0; noConnections=0;
			break;
	}

   if(openptr==NULL){
      fprintf(stderr,"The file did not open\n");
	  dfl2=1; /*show TRP-Cage as default*/
      return;
   }

	fgets(firstField,7,openptr);

   /*go til firstField is “ATOM”*/
	while(strcmp(firstField,"ATOM  ")!=0 && strcmp(firstField,"HETATM")!=0 && kk<100000){
		fgets(blank,1024,openptr);
		fgets(firstField,7,openptr);kk++;
	}/*Now firstField == ATOM*/
dummyNums[0][1]=kk;

   while( strcmp(firstField,"END   ")!=0 && strcmp(firstField,"ENDMDL")!=0 && kk<100000){
      chaiNum++;	jj=1;	kk++; dummyNums[0][2]=kk; strcpy(dummyChars[0][0],firstField);

	while( strcmp(firstField,"TER   ")!=0 && kk<100000){
			kk++; dummyNums[0][3]=kk; strcpy(dummyChars[0][1],firstField);


       	if(strcmp(firstField,"ATOM  ")==0 || strcmp(firstField,"HETATM")==0){
			ii++; /*line number*/ strcpy(dummyChars[0][2],firstField);


	/*use ATOM/HETATM record format*****/
		fgets(serialNum,6, openptr);/*7-11*/
		atomInfoNum[ii][ATONUM]=atoi(serialNum);
		fgets(blank,2,openptr);/*12*/
		fgets(atomName,5,openptr);/*13-16*/
		strcpy(atomInfoStr[ii][CONINF],atomName);
		fgets(altLoc,2,openptr);/*17*/
		fgets(resName,4,openptr);/*18-20*/
		strcpy(atomInfoStr[ii][RESTYP],resName);
		fgets(blank,2,openptr);/*21*/
		fgets(chainID,2,openptr);/*22*/
		strcpy(atomInfoStr[ii][CHANID],chainID);
		fgets(resNum,5,openptr);/*23-26*/
		atomInfoNum[ii][RESNUM]=atoi(resNum);
		fgets(iCode,2,openptr);/*27*/
		fgets(blank,4,openptr);/*28-30*/
		fgets(coord,9,openptr);/*31-38*/
		aCent[ii][X]=atof(coord);
		fgets(coord,9,openptr);/*39-46*/
		aCent[ii][Y]=atof(coord);
		fgets(coord,9,openptr);/*47-54*/
		aCent[ii][Z]=atof(coord);
		fgets(occup,7,openptr);/*55-60*/
		fgets(tempf,7,openptr);/*61-66*/
		fgets(blank,7,openptr);/*67-72*/
		fgets(segID,5,openptr);/*73-76*/
		fgets(elementName,3,openptr);/*77-78*/
		strcpy(atomInfoStr[ii][ATOTYP],elementName);
		/*************fgets(charge,3,openptr);/*79-80*/
		fgets(blank,1024,openptr);/*rest of line*/

			fgets(firstField,7,openptr);


			/********************
			if(virg){virg=0;strcpy(lastChainID,chainID);}
			/************
			if(strcmp(lastChainID,chainID)!=0){nnnn++;strcpy(lastChainID,chainID);}
			
			  /************************
			if(strcmp(lastChainID,chainID)!=0){
				nnnn++;
				chaiNum++;
				begChainInx[chaiNum]=ii;
				strcpy(lastChainID,chainID);
			}
			/*************************************/

		/*translate atom type to number*/
			if(strcmp(elementName," H")==0){atomInfoNum[ii][ATOTYP]=HYDR;}
			else{
			if(strcmp(elementName," C")==0){atomInfoNum[ii][ATOTYP]=CARB;}
			else{
       		if(strcmp(elementName," N")==0){atomInfoNum[ii][ATOTYP]=NITR;}
			else{
			if(strcmp(elementName," O")==0){atomInfoNum[ii][ATOTYP]=OXYG;}
			else{
			if(strcmp(elementName," S")==0){atomInfoNum[ii][ATOTYP]=SULP;}
			else{
			atomInfoNum[ii][ATOTYP]=NULL;
			}}}}}
		/*record indices for backbone and alpha carbon atoms */
       		if(strcmp(atomName," N  ")==0){bbIndices[chaiNum][jj]=ii;jj++;}
			if(strcmp(atomName," CA ")==0){bbIndices[chaiNum][jj]=ii;jj++;
       			alphaCIndices[aCInx]=ii;lastResidue++; aCInx++;}
			if(strcmp(atomName," C  ")==0){bbIndices[chaiNum][jj]=ii;jj++;}

		/*translate residue type to number*/
       if(strcmp(resName,"GLY")==0){atomInfoNum[ii][RESTYN]=GLY;}
	else{
	if(strcmp(resName,"ALA")==0){atomInfoNum[ii][RESTYN]=ALA;}
	else{
	if(strcmp(resName,"VAL")==0){atomInfoNum[ii][RESTYN]=VAL;}
	else{
	if(strcmp(resName,"LEU")==0){atomInfoNum[ii][RESTYN]=LEU;}
       else{if(strcmp(resName,"ILE")==0){atomInfoNum[ii][RESTYN]=ILE;}
	else{
	if(strcmp(resName,"MET")==0){atomInfoNum[ii][RESTYN]=MET;}
       else{
      if(strcmp(resName,"PHE")==0){atomInfoNum[ii][RESTYN]=PHE;}
	else{
      if(strcmp(resName,"PRO")==0){atomInfoNum[ii][RESTYN]=PRO;}
	else{
      if(strcmp(resName,"TRP")==0){atomInfoNum[ii][RESTYN]=TRP;}
	else{
      if(strcmp(resName,"ARG")==0){atomInfoNum[ii][RESTYN]=ARG;}
	else{
      if(strcmp(resName,"ASP")==0){atomInfoNum[ii][RESTYN]=ASP;}
	else{
      if(strcmp(resName,"GLU")==0){atomInfoNum[ii][RESTYN]=GLU;}
	else{
      if(strcmp(resName,"LYS")==0){atomInfoNum[ii][RESTYN]=LYS;}
	else{
      if(strcmp(resName,"HIS")==0){atomInfoNum[ii][RESTYN]=HIS;}
       else{
      if(strcmp(resName,"ASN")==0){atomInfoNum[ii][RESTYN]=ASN;}
	else{
      if(strcmp(resName,"CYS")==0){atomInfoNum[ii][RESTYN]=CYS;}
	else{
      if(strcmp(resName,"GLN")==0){atomInfoNum[ii][RESTYN]=GLN;}
	else{
       if(strcmp(resName,"SER")==0){atomInfoNum[ii][RESTYN]=SER;}
	else{
      if(strcmp(resName,"THR")==0){atomInfoNum[ii][RESTYN]=THR;}
	else{
      if(strcmp(resName,"TYR")==0){atomInfoNum[ii][RESTYN]=TYR;}
       }}}}}}}}}}}}}}}}}}}/*end 19th else: translating resName(char) to restyn(int)*/

		}/*end if for ATOM/HETATM line*/
		else{/*something is seriously wrong/flush the line*/
			fgets(blank,1024,openptr);
			fscanf(openptr,"%s",firstField);kk++; dummyNums[0][4]=kk;
		}

	}/*end while(!=TER): we have reached a TER; increment chaiNum & flush line*/
	/*virg=1;/*qwert:reset for next pass through wile loop*/
	begChainInx[chaiNum+1]=ii+1;
	fgets(blank,1024,openptr);
	fscanf(openptr,"%s",firstField);
	lastBB[chaiNum]=jj-1;

	/**While we are not at END, ENDMDL, or ATOM/HETATM, flush**/
	while(strcmp(firstField,"ATOM  ")!=0 && strcmp(firstField,"HETATM")!=0 && strcmp(firstField,"END   ")!=0 && strcmp(firstField,"ENDMDL")!=0 && kk<100000){
		fgets(blank,1024,openptr);
		fgets(firstField,7,openptr);kk++; dummyNums[0][3]=kk;
strcpy(dummyChars[0][3],firstField);
dummyNums[0][9]=strcmp(firstField, "END   ");
dummyNums[0][10]=strcmp(firstField, "END  ");
dummyNums[0][11]=strcmp(firstField, "END ");
dummyNums[0][12]=strcmp(firstField, "END");
	}

	}/* end while(!=END/ENDMDL)*/


	begChainInx[1]=1;
	lastAtom=ii;
	lastChain=chaiNum;
	centerTheModel();
	computeBackBone();
	/**computeFrame();**/
	/*constructMolecule();*/
	reConstructMarker();
}/*end openpdbfile()*/

/**********************************************************************/
/**********************************************************************/
/* Stuart Levy  1998  improved gadgets  getnumber, bump               */
/**********************************************************************/
float getnumber(float dflt){    /* dflt is an action or a constant */
      if(!hasnumber)return dflt;  /* in case no number to put in */
      tmp = sign ? -number : number;
      return decimal>0 ? tmp/(float)decimal : tmp ;
}
/**********************************************************************/
void bump(float *val, float incr){ /* use bump(&nose, .01) for example */
     float abs = fabs(incr); char  frmt[8], code[32]; int digits = 1;
     if(hasnumber){ *val = getnumber(0); return; }
     if(abs <= .003) digits=3; else if(abs <= .03) digits=2;
     sprintf(frmt, "%%.%de" , digits);
     sprintf(code, frmt, incr >0 ? (*val)*(1+abs):(*val)/(1+abs));
     sscanf(code, "%f", val);
}
/**********************************************************************/
/* keyboard - handles input from the keyboard.                       */
/*   I needed more keys, so I commented out some of the classic skel  */
/*   state variables and stole their keys.  ~baf                      */
/**********************************************************************/
void keyboard(unsigned char key, int x, int y){
     #define  IF(K)            if(key==K)
     #define  PRESS(K,A,b)     IF(K){b;}IF((K-32)){A;}
     /*was backwards in previous versions */
     #define  TOGGLE(K,flg)    IF(K){(flg) = (int)getnumber(1-(flg));}
     #define  CYCLE(K,f,m)     PRESS((K), \
              getnumber((f)=(((f)+(m)-1)%(m))), getnumber((f)=(++(f)%(m))))
   /* Only ASCII characters can be processes by this GLUT callback function */
   IF(27)exit(0);                                         /* ESC exit            */
   CYCLE(' ', mode,FLYMODE+1);                           /* fly/turn modes      */
   CYCLE('w',msg,3);                                      /* writing on/off      */
   PRESS('s', bump(&speed,-.02), bump(&speed,.02))        /* flying speed        */
   PRESS('q', bump(&torq, -.02),  bump(&torq,.02))        /* turning speed*/
   PRESS('z', dfl2=1;/*deFault();reSet()*/, reSet())                 /* zap changes         */
   TOGGLE('p', friz)									  /* stop spinning */
   TOGGLE('x', axx)
   PRESS('y', wfar *= 1.01,   wfar   /= 1.01) /* rear clipping plane*/
   PRESS('u', focal *= 1.01,   focal   /= 1.01) /* rear clipping plane*/

	/**development elements**/
    TOGGLE('h', whiteBG)
    CYCLE('/', chapmat,4)
	/**PRESS('/', checkProcess(), checkProcess())/**/
   	PRESS('?', recordAFF(), recordAFF())
   	/**end development elements**/
	if(sceneType==MOLECULE){
   		PRESS('a', bump(&fac1,-.02), bump(&fac1,.02)) /* multiplier for forward vector */
		CYCLE('b', backBone, LASTBBONE+1)
  		PRESS('c', if(targetAC>1) targetAC--; reConstructMarker(), if(targetAC<lastResidue) targetAC++; reConstructMarker());   /* atom to mark */
   		PRESS('f', bump(&SF,-.2);setAtomRadii(), bump(&SF,.2);setAtomRadii())        /*  */
 /*  		CYCLE('i', openFileNum,8)
 */
		PRESS('i',changeMol=1,changeMol=1);
   		CYCLE('j', whichChain,20)
   		TOGGLE('l', selectChains)
   		PRESS('m', TOGGLE('m', bigMarker);reConstructMarker(), TOGGLE('m', bigMarker);reConstructMarker());   /* change size of marker */
		TOGGLE('n', onFly)
/*   		PRESS('o', openpdbfile(),openpdbfile());   /* atom to mark */
   		TOGGLE('r', fram1)
	}


   	PRESS('d', if(cMax>2) cMax /= 2, cMax *= 2);   /* animation speed */
   	TOGGLE('e', changeCent)
 /***
 PRESS('g',TOGGLE('g',geo);if(geo)binoc=-1;else binoc=0;,TOGGLE('g',geo);if(geo)binoc=-1;else binoc=0;)
 */
 PRESS('g',CYCLE('g',geo,3);if(geo==0){binoc=0;} if(geo==1){binoc=1;}if(geo==2){binoc=-1;},CYCLE('g',geo,3);if(geo==0){binoc=0;} if(geo==1){binoc=1;}if(geo==2){binoc=-1;})

	PRESS('n', if(nn1>0) nn1--,  nn1++);   /* atom to mark */
 /*  	PRESS('t',CYCLE('t', sceneType,5);reSet(),CYCLE('t', sceneType,5);reSet());  */
	PRESS('t',changeScene=1,changeScene=1);
 /*  	TOGGLE('u', opaque)
 */
	TOGGLE('v', mov)

   /* Stuart Levy's gadget parser from avn.c */
   if(key >= '0' && key <= '9'){ hasnumber = 1;
         number = 10*number + (key - '0'); decimal *= 10;}
   else if(key == '.'){ decimal = 1;}
   else if(key == '-'){ sign = -1;}
   else {hasnumber = number = decimal = sign = 0;}

#ifndef SZG
   glutPostRedisplay();
#endif
}
/**********************************************************************/
/* special_keybo - handles non-ASCII keypresses                       */
/**********************************************************************/
void special_keybo(int key, int x, int y){
/*non-ASCII keypresses go here, if you're lucky enough to know their names */
 fprintf(stderr," non-ASCII character was pressed.\n");
 fprintf(stderr," use special_keybo() to process it\n");
}
/**********************************************************************/
/* speedometer - Calculates framerate                                 */
/**********************************************************************/
#ifdef WIN32
float speedometer(void){                      /* this one is for win32*/
  double dbl; static double rate; static int ii=0;
  static struct _timeb lnow, lthen;
    if(++ii % 8 == 0){                 /* 8 times around measure time */
	   _ftime(&lnow);
		dbl =  (double)(lnow.time - lthen.time)
			 +(double)(lnow.millitm - lthen.millitm)/1000;
		lthen = lnow;  rate = 8/dbl;
      }
    return((float)rate);
}
#endif
/**********************************************************************/
#ifdef IRIX
float speedometer(void){
double dbl; static double rate; static int ii=0;
static struct timezone notused; static struct timeval now, then;
   if(++ii % 8 == 0){  /* 8 times around measure time */
      gettimeofday(&now, &notused); /* elapsed time */
      dbl =  (double)(now.tv_sec - then.tv_sec)
         +(double)(now.tv_usec - then.tv_usec)/1000000;
      then = now;  rate = 8/dbl;
      }
   return((float)rate);
}
#endif
//#ifdef OSXMac
#ifdef LINUX
float speedometer(void){
double dbl; static double rate; static int ii=0;
static struct timezone notused; static struct timeval now, then;
   if(++ii % 8 == 0){  /* 8 times around measure time */
      gettimeofday(&now, &notused); /* elapsed time */
      dbl =  (double)(now.tv_sec - then.tv_sec)
         +(double)(now.tv_usec - then.tv_usec)/1000000;
      then = now;  rate = 8/dbl;
      }
   return((float)rate);
}
#endif


/**********************************************************************/
/* char2wall - ??                                                     */
/**********************************************************************/
void char2wall(float x,float y,float z, char buf[]){
     char *p; glRasterPos3f(x,y,z);
     for(p = buf;*p;p++) glutBitmapCharacter(GLUT_BITMAP_9_BY_15,*p);
}
/**********************************************************************/
/* messages - writes text on the screen                               */
/**********************************************************************/
void messages2(void){
 char buf[256];
/**********/
#define  LABEL2(x,y,W,u) {sprintf(buf,(W),(u));char2wall(x,y,0.,buf);}
/************/
char phrase[256];
#define  LABEL(x,y,W,u) sprintf(phrase,(W),(u)); cmov2((x),(y)); charstr(phrase)

  glMatrixMode(GL_PROJECTION); glPushMatrix(); glLoadIdentity();
  gluOrtho2D(0,3000,0,3000);
  glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();


         if(mode==TURNMODE) glColor3f(0x22/255.,0x88/255.,0xdd/255.);
         else      glColor3f(.8,.8,.8);

/***********************
if(geo){LABEL2(750,1500,"%s","o");}

	else{LABEL2(1499,1500,"%s","o");}
/*************************/

	/** signal turning mode by text color **/
		 if(mode==TURNMODE) glColor3f(1, 0, 1);
         else      glColor3f(1, 1, 0);

LABEL2(80,80,"%4.1f fps",speedometer());
if(virg){
	LABEL2(80,2840,"Welcome! Press the w-key several times. %s","");
	LABEL2(500,10,"Navigating in 3D Models %s","");
	if(msg!=1){virg=0;}
}
else{/*the usual messages, following “Welcome”*/
	switch(msg){
		case VERBOSE:/*==2**/
 	  		LABEL2(80,2500,"(S)peed  %0.3f",speed);
      		LABEL2(80,2400,"tor(Q) %0.3f",torq);
      		if(sceneType!= MOLECULE){LABEL2(80,2300,"animSp(d) %d",cMax);}
			else{LABEL2(80,2300,"Molecule type: %s", showMolName); }
      		LABEL2(80,2200,"(Esc)ape %s","");
      		LABEL2(80,2100,"(p)ause %s","");
      		LABEL2(80,2000,"(g)eowall %s","");
      		LABEL2(80,1900,"(Z)ap %s","");
 	     	LABEL2(80,1800,"mo(v) %d",mov);
/* 	     	LABEL2(80,1700,"(/)chapmat %d",chapmat);
*/
LABEL2(80,1700,"opeFileNum %d",openFileNum);
LABEL2(80,1600,"nnnn %d",nnnn);
/*qwer*/

      		if(sceneType== MOLECULE){
      			LABEL2(80,900,"chgC(e)nt %d",changeCent);
				switch(openFileNum){
					case 0:LABEL2(80,800,"wh(i)ch file: %s","0. E. coli subunit C");break;
					case 1:LABEL2(80,800,"wh(i)ch file: %s","1. E. coli N-Term Dom. delta");break;
					case 2:LABEL2(80,800,"wh(i)ch file: %s","2. E. coli subunit B");break;
					case 4:LABEL2(80,800,"wh(i)ch file: %s","3.  spinach F1");break;
					case 3:LABEL2(80,800,"wh(i)ch file: %s","4.  yeast rotor motor ");break;
					case 5:LABEL2(80,800,"wh(i)ch file: %s"," CRP");break;
					case 6:LABEL2(80,800,"wh(i)ch file: %s"," CRP from VMD");break;
					case 7:LABEL2(80,800,"wh(i)ch file: %s"," x.pdb");break;
				}
      			LABEL2(80,700,"(o)pen file %s","");
      			LABEL2(80,600,"(m)arker size %s","");
      			LABEL2(80,500,"(f)ill atom %0.3f",SF);
      			LABEL2(80,400,"se(l)ectChain %d",selectChains);
      			LABEL2(80,300,"f(r)ame %d",fram1);
     			if(selectChains){LABEL2(500,400,"(j) %d",whichChain);}
			}
    			LABEL2(10,10,"30jul06 11jun07 %s","");
    			LABEL2(1280,10,", P McCreary, G Francis, D Bergman, UIUC+Evergreen %s","");
			/** no break, just plow ahead to catch the always stuff **/

		case SPARCE:
		    	LABEL2(500,10,"zatp %s","");
 			LABEL2(80,2840,\
	      		"(W)riting (MAUS)Fore (BAR)Flymode: %s",
             	mode?"FLYING":"CONTROL");

			switch(sceneType){
				case LOXODROMIC:
		   			LABEL2(80,2700,"(t)ype: Infinite Cyclic, Loxodromic",0);
		 			break;
				case ELLIPTIC:
		    			LABEL2(80,2700,"(t)ype: Finite Cyclic, Elliptic",0);
		 			break;
				case KLEIN4:
		    			LABEL2(80,2700,"(t)ype: Klein Four",0);
					break;
				case PARABOLIC:
		    			LABEL2(80,2700,"(t)ype: Parabolic",0);
		 			break;
				case MOLECULE:
		    			LABEL2(80,2700,"(t)ype: Molecule Model; (b)ackbone %i",backBone);
						LABEL2(800,2600,"amino a(c)id #%i",targetAC);
      					LABEL2(80,2600,"residue type %s", atomInfoStr[alphaCIndices[targetAC]][RESTYP]);
					break;
			}
	}
}/**end else for regular messages following “Welcome”**/

	glPopMatrix();
	glMatrixMode(GL_PROJECTION); glPopMatrix();

}
void messagesOld(void){
  char buf[256]; /* console messages are done differently from cave */
#define  LABEL2(x,y,W,u) {sprintf(buf,(W),(u));char2wall(x,y,0.,buf);}
  glMatrixMode(GL_PROJECTION); glPushMatrix(); glLoadIdentity();
  gluOrtho2D(0,3000,0,3000);
  glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();

    /*bull's eye*/
         if(mode==TURNMODE) glColor3f(0x22/255.,0x88/255.,0xdd/255.);
         else      glColor3f(.8,.8,.8);

	if(geo){LABEL2(750,1500,"%s","o");}

	else{LABEL2(1499,1500,"%s","o");}



	  /* writings */
         if(mode==TURNMODE) glColor3f(1, 0, 1);
         else      glColor3f(1, 1, 0);


      LABEL2(80,80,"%4.1f fps",speedometer());
      LABEL2(80,2840,\
      "(ESC)ape (MAUS2)Fore (BAR)Flymode %s (W)riting (P)ause",
             mode?"FLYING":"CONTROL");
      LABEL2(10,10,"8 June 2004. Elementary Kleinian Groups.  Paul R. McCreary, Xavier University of LA %s","");
      LABEL2(80,2700,"(S)peed  %0.3f",speed);
      LABEL2(80,2650," tor(Q) %0.3f",torq);
      LABEL2(80,2600,"(Z)ap %s","");
      LABEL2(80,2550,"(T)ime = %1i", tt);

	if(geo){LABEL2(80,2500, "(g)eowall %i", geo);}
	if(!geo){LABEL2(80,2500, "bi(N)oc %i", binoc);}

	  	switch(sceneType){
		case LOXODROMIC:
   LABEL2(80,2400,"Grp(T)ype=Infinite Cyclic, Loxodromic",0);
   LABEL2(80,2350,"# it(e)rates=%d",numIterates);
   LABEL2(80,2300,"P(a)th",0);
   LABEL2(80,2250,"Tw(i)st =Pi*%g",twist/Pi);
   LABEL2(80,2200,"(K)multiplier = %g",multiplierK);
   LABEL2(80,2150,"(L)ayers",0);
		break;
		case ELLIPTIC:
    LABEL2(80,2400,"Grp(T)ype=Finite Cyclic, Elliptic",0);
    LABEL2(80,2350,"aniSp(d) %d",cMax);
    LABEL2(80,2300,"(O)rder=%d",order);
    LABEL2(80,2250,"(M)ov",0);
    LABEL2(80,2200,"(H)oro=%d",horo);
		break;
		case KLEIN4:
    LABEL2(80,2400,"Grp(T)ype=Klein Four",0);
    LABEL2(80,2350,"s(c)ene=%d",sceneNum);
    LABEL2(80,2300,"aniSp(d) %d",cMax);
    LABEL2(80,2250,"(M)ov",0);

		break;
		/**case ELLOX:
    LABEL2(80,2400,"Grp(T)ype=Loxodromic X Elliptic",0);
    LABEL2(80,2350,"(F1)# lox appl=%d",loxCount);
    LABEL2(80,2300,"(F5)# ell appl=%d",ellCount);
    LABEL2(80,2250,"aniSp(d) %d",cMax);
    LABEL2(80,2200,"c(h)ase",0);
    LABEL2(80,2150,"thisOne=%d",thisOne);

		break;**/
		case PARABOLIC:
    LABEL2(80,2400,"Grp(T)ype=Parabolic",0);
    LABEL2(80,2350,"(arrows) slidePnt=%d",slideInx);
		break;
	}/**qaz**/

	LABEL2(80,370,"mov=%d",mov);
	LABEL2(80,320,"loxLength=%d",loxLength);
	LABEL2(80,270,"lMax=%d",lMax);
   LABEL2(80,220,"movieCounter=%d",movieCounter);

    LABEL2(80,170,"P(l)ane (B)ands",0);

	  glPopMatrix();
glMatrixMode(GL_PROJECTION); glPopMatrix();

}
/**********************************************************************/
/* chaptrack - Handles mouse input durring FLYMODE/TURNMODE           */
/**********************************************************************/
void chaptrack(int but,int xx,int yy,int shif){
   long dx,dy;

	   if(geo){dx = xx -.25*xt; dx = abs(dx)>5?dx:0;}
	else{ dx = xx -.5*xt; dx = abs(dx)>5?dx:0;}

   dy = yy -.5*yt; dy = abs(dy)>5?dy:0;

   if(friz /**|| showslate**/) return; /* freeze if we are in slate mode or friz */

   glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();
   pushit(trans);
   if(mode==TURNMODE){glTranslatef(trans[X],trans[Y],trans[Z]);} /**//****/
   glRotatef(dx*torq,0.,1.,0.); glRotatef(dy*torq,1.,0.,0.);
   if(mode==FLYMODE){
      glPushMatrix();
      glMultMatrixf(starmat);
      glGetFloatv(GL_MODELVIEW_MATRIX,starmat);
      glPopMatrix();
   }
#ifndef OSXMac
   if(but&(1<<GLUT_MIDDLE_BUTTON))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(but&(1<<GLUT_RIGHT_BUTTON ))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(but&(1<<GLUT_LEFT_BUTTON  ))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(flyForward)glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
#else
   if(but&(1<<GLUT_LEFT_BUTTON)) glTranslatef(0.,0.,shif ? -speed : speed);
#endif
  if(mode==TURNMODE){glTranslatef(-trans[X],-trans[Y],-trans[Z]);} /**//****/
   glMultMatrixf(aff);
   glGetFloatv(GL_MODELVIEW_MATRIX,aff);
   FOR(ii,0,3){luxx[ii]=0; FOR(jj,0,3)luxx[ii] +=aff[ii*4+jj]*lux[jj];}
   glPopMatrix();
}
void chaptrack2(int but,int xx,int yy,int shif){
   long dx,dy;

	   if(geo){dx = xx -.25*xt; dx = abs(dx)>5?dx:0;}
	else{ dx = xx -.5*xt; dx = abs(dx)>5?dx:0;}

   dy = yy -.5*yt; dy = abs(dy)>5?dy:0;

   if(friz /**|| showslate**/) return; /* freeze if we are in slate mode or friz */

   glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();
   pushit2(trans);
   if(mode==TURNMODE){glTranslatef(trans[X],trans[Y],trans[Z]);} /**//****/
   glRotatef(dx*torq,0.,1.,0.); glRotatef(dy*torq,1.,0.,0.);
   if(mode==FLYMODE){
      glPushMatrix();
      glMultMatrixf(starmat);
      glGetFloatv(GL_MODELVIEW_MATRIX,starmat);
      glPopMatrix();
   }
#ifndef OSXMac
   if(but&(1<<GLUT_MIDDLE_BUTTON))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(but&(1<<GLUT_RIGHT_BUTTON ))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(but&(1<<GLUT_LEFT_BUTTON  ))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(flyForward)glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
#else
   if(but&(1<<GLUT_LEFT_BUTTON)) glTranslatef(0.,0.,shif ? -speed : speed);
#endif
  if(mode==TURNMODE){glTranslatef(-trans[X],-trans[Y],-trans[Z]);} /**//****/
   glMultMatrixf(aff2);
   glGetFloatv(GL_MODELVIEW_MATRIX,aff2);
   FOR(ii,0,3){luxx[ii]=0; FOR(jj,0,3)luxx[ii] +=aff2[ii*4+jj]*lux[jj];}
   glPopMatrix();
}/*end chaptrack2*/
void chaptrack3(int but,int xx,int yy,int shif){
   long dx,dy;

	   if(geo){dx = (long)(xx -.25*xt); dx = (long)(abs(dx)>5?dx:0);}
	else{ dx = (long)(xx -.5*xt); dx = (long)(abs(dx)>5?dx:0);}

   dy = (long)(yy -.5*yt); dy = (long)(abs(dy)>5?dy:0);

   if(friz /**|| showslate**/) return; /* freeze if we are in slate mode or friz */


   glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();
   pushit3(trans);
   if(mode==TURNMODE){glTranslatef(trans[X],trans[Y],trans[Z]);} /**//****/
   glRotatef(dx*torq,0.,1.,0.); glRotatef(dy*torq,1.,0.,0.);
   if(mode==FLYMODE){
      glPushMatrix();
      glMultMatrixf(starmat);
      glGetFloatv(GL_MODELVIEW_MATRIX,starmat);
      glPopMatrix();
   }
#ifndef OSXMac
   if(but&(1<<GLUT_MIDDLE_BUTTON))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(but&(1<<GLUT_RIGHT_BUTTON ))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(but&(1<<GLUT_LEFT_BUTTON  ))glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
   if(flyForward)glTranslatef(0.,0.,shif ? -speed*speedMed : speed*speedMed);
#else
   if(but&(1<<GLUT_LEFT_BUTTON)) glTranslatef(0.,0.,shif ? -speed : speed);
#endif
  if(mode==TURNMODE){glTranslatef(-trans[X],-trans[Y],-trans[Z]);} /**//****/
   glMultMatrixf(aff2);
   glGetFloatv(GL_MODELVIEW_MATRIX,aff2);
   FOR(ii,0,3){luxx[ii]=0; FOR(jj,0,3)luxx[ii] +=aff2[ii*4+jj]*lux[jj];}
   glPopMatrix();
}/*end Chaptrack3*/

/**********************************************************************/
/* viewfrom - Handles the changes for riding in the cockpit           */
/**********************************************************************/
void viewfrom(int lcockpit){  /* was just "cockpit" but CAVE macro failed here */
    if(!lcockpit){
#ifndef OSXMac
      // drawstars(starmat);
#endif
      if(binoc) glTranslatef(-binoc*nose,0.0,0.0);
      /******
      switch(chapmat){
		  case 0: glMultMatrixf(aff);break;
		  case 1: glMultMatrixf(aff2);break;
	  }/********/
   }/*end if(!1cockpit)*/
   else{
    float inv[16];
                       /* if there is showchaseplane, view from it */
             /**invert(inv, showchaseplane? fram5: fram3); **/

       glTranslatef(0.,-BIGCAR,0.); /* moves the head above the track */
      glMultMatrixf(inv);
      glPushMatrix();
         glScalef(9,9,9);  /*makes sure stars stay a good size*/
#ifndef OSXMac
         // drawstars(starmat);
#endif
      glPopMatrix();
   }/*end else for if(!1cockpit)*/

}
/**********************************************************************/
/**********************************************************************/
/* reshaped - handles window resizing                                 */
/**********************************************************************/
void reshaped(int xx, int yy){xt=xx ; yt=yy;}   /* origin of moved window */
/**********************************************************************/
/* drawcons - draws everything                                        */
/*   - Now uses showslate variable to draw multiple viewports         */
/**********************************************************************/
void drawcons(void){
   if(whiteBG){glClearColor(0.950f,0.950f,0.950f,1.0f);}
   glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
   dawn(1);

    if(binoc){

			if(geo){glViewport(0,0,xt/2,yt);}

			else{glViewport(0,yt/4,xt/2,yt/2);}

	}
     else  glViewport(0,0,xt,yt);

      /***********Main*Window*******************/
      glMatrixMode(GL_PROJECTION); glLoadIdentity();

      	if(geo){glFrustum(-mysiz*xt/yt,mysiz*xt/yt,-2*mysiz,2*mysiz,mysiz*focal,wfar);}

		else{glFrustum(-mysiz*xt/yt,mysiz*xt/yt,-mysiz,mysiz,mysiz*focal,wfar); }

      glMatrixMode(GL_MODELVIEW); glLoadIdentity();
      viewfrom(cockpit);
      /**drawall();/**/drawallllll();
      if(binoc){

			if(geo){glViewport(xt/2,0,xt/2,yt);}

			else{glViewport(xt/2,yt/4,xt/2,yt/2);}

         glMatrixMode(GL_MODELVIEW);
         glLoadIdentity();
#ifndef OSXMac
         // drawstars(starmat);
#endif
         glTranslatef(binoc*nose,0.0,0.0);
         switch(chapmat){/*?????????????????????*/
			 case 0: glMultMatrixf(aff);
			 		drawall(); break;
			 case 1: glMultMatrixf(aff2);
			 		drawall2(); break;
		 }
         /**drawall();/**/
      }/*end if(binoc)*/

   glViewport(0,0,xt,yt);
   if(msg) messages2();if(mov) movie();
   glutSwapBuffers();

}
/**********************************************************************/
/* idle - function that gets executed when the CPU is idle            */
/**********************************************************************/
void idle(void){ /*do this when nothing else is happening*/
   if(morph) autotymer(0);  /* advance autotymer */
   glutPostRedisplay();  /*redraw the window*/
   switch(chapmat){
	   case 0: chaptrack(BUT,XX,YY,SHIF); break;
	   case 1: chaptrack2(BUT,XX,YY,SHIF); break;
	   case 2:  chaptrack(BUT,XX,YY,SHIF);  chaptrack2(BUT,XX,YY,SHIF); break;
	   case 3: chaptrack(BUT,XX,YY,SHIF); chaptrack3(BUT,XX,YY,SHIF);break;
   }
   glMultMatrixf(aff);findscreen();

	/* operations carried out in postexchange, and so should be done for console also */
    if(dfl2){dfl2=0;deFault();reSet();}
	if(changeScene){
		changeScene=0;
		if(sceneType<4){sceneType=sceneType+1;} 
		else {sceneType=0;};
		reSet();
	}
	if(changeMol){
		changeMol=0;
		if(openFileNum<(maxMol-1)){openFileNum=openFileNum+1;} 
		else {openFileNum=0;};
		openpdbfile();
	}/*** End console copy of postexchange operations for console****/
}
/**********************************************************************/
/* mousepushed - mouse button handler. I use this to determine dragmode*/
/**********************************************************************/
void mousepushed(int but,int stat,int x,int y){
  if(stat==GLUT_DOWN) BUT |= (1<<but);
  else BUT &= (-1 ^ (1<<but));
  XX=x; YY=y; SHIF=(glutGetModifiers()==GLUT_ACTIVE_SHIFT)?1:0;

}
/**********************************************************************/
/* mousemoved - mouse movement handler                                */
/**********************************************************************/
void mousemoved(int x,int y){  XX=x; YY=y; }
/**********************************************************************/
/* dataprep - intializes data                                         */
/**********************************************************************/
void dataprep(void){
 	deFault();
	reSet();
#ifndef OSXMac
        initstars();
#endif
        
}
/**********************************************************************/
/* %%%%%%%%%%%%%      CAVEISH        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/**********************************************************************/

#ifdef CAVE
void cavekeybo(void){
#define  CAVEIF(K)         if(CAVEgetbutton(K))
#define  CAVESOAK(K)        while(CAVEgetbutton(K))

#define  CAVEPRESS(K,A,b)  CAVEIF(K){if(CAVEgetbutton(CAVE_LEFTSHIFTKEY)||\
                           CAVEgetbutton(CAVE_RIGHTSHIFTKEY)){b;}else{A;}\
                           CAVESOAK(K);}
#define  CAVETOGGLE(K,flg) CAVEIF(K){(flg) = 1-(flg); CAVESOAK(K);}
#define  CAVECYCLE(K,f,m)  CAVEPRESS((K),(f)=(((f)+(m)-1)%(m)),(f)=(++(f)%(m)) )
#define  CAVESLIDI(K,f,m,M)  CAVEPRESS(K,(--f<m?m:f), (++f>M?M:f))
#define  CAVESLIDF(K,f,m,M,d) CAVEPRESS(K,((f -= d)<m?m:f), ((f += d)>M?M:f))
CAVEPRESS(CAVE_PKEY, wfar *= 1.01,   wfar   /= 1.01) /* rear clipping plane*/
CAVEPRESS(CAVE_QKEY, speed /= 1.02, speed *= 1.02); /* flying speed       */
CAVEPRESS(CAVE_RKEY, torq /= 10.02,  torq *= 10.02); /* turning speed      */
CAVETOGGLE(CAVE_YKEY,wnd);  /* mauspaw vs wandpaw *//* not yet implemented */
CAVEPRESS(CAVE_ZKEY, deFault(), deFault());        /* zap changes       */
CAVETOGGLE(CAVE_FKEY, showframe );
CAVETOGGLE(CAVE_UKEY, showcube );
CAVETOGGLE(CAVE_KKEY,  cockpit  );


CAVETOGGLE(CAVE_DKEY, showbishop );
CAVETOGGLE(CAVE_YKEY, autopilot );
CAVETOGGLE(CAVE_EKEY, showchaseplane );
} /* end cavekeybo */

void graffiti(void){             /* used to be called speedo  */
  char buf[256];                 /*messages written into here */
#define  LABEL3(x,y,z,W,u){sprintf(buf,(W),(u));char2wall(x,y,z,buf);}
  static float last=0.,fps;      /*for measuring time         */
  if(!CAVEMasterDisplay()) return;  /* USE CAVEMasterWall in papeBeta */
                                    /* in grafiXlab we have old CAVElib? */
  if(floor(2.*(*CAVETime))>floor(2.*last)){  /* rationalize this */
     last=*CAVETime;
     fps=*CAVEFramesPerSecond;
     }
  if(mode==TURNMODE){
        glColor3f(1.,0.,1.);
        LABEL3(-3.8,1.,-5.0,\
        "4=auto 2=mode 3=shrink 6=grow 5=freeze 7=reset %s","");
        } else {
        glColor3f(1.,1.,0.);
        LABEL3(-3.8,1.0,-5.0,
          " push thumb-button to fly, 2=tractor gap %s","gap");
        }
     LABEL3(-4.8,8.0,-5.0, "bishopCoaster\
    by Ben Farmer and Karan Dalal, im2001 (C)Summer 2001  U.Illinois %s","")
       LABEL3(-4.8,1.0,-5.0,"%5.1f fps",fps);
       LABEL3(-4.8,7.5,-5.0,"Time = %i",tt);
       LABEL3(-4.8,7.0,-5.0,"d(T) = %3i",dtt);
}

void cavetrack(void){
  float azi,ele,rol,hx,hy,hz,wx,wy,wz;
  static float ohx,ohy,ohz;   /*old head position */
  static float owx,owy,owz;   /*old wand position */
  static int opaw = 0,paw, joy = 0, friz = 0;

  if(!CAVEMasterDisplay()) return;  /* doit once per frame */
  paw = wnd ? (CAVEBUTTON1*2+CAVEBUTTON2)*2+CAVEBUTTON3 : mauspaw;
  if(paw) joy = 1; /* activate joystick by clicking any button*/
  if(opaw!=2 && paw==2) mode = (mode+1)%(TURNMODE+1); /* click modes */
  if(mode==TURNMODE){   /* hold middle button and click right to shrink */
     if (opaw == 2 && paw == 3)siz /= 1.5;     /*shrink object */
     if (opaw == 2 && paw == 6)siz *= 1.5;     /* grow   object */
     if (opaw != 4 && paw == 4)cockpit=1-cockpit; /* viewchanger */
     if (opaw != 5 && paw == 5)friz = 1-friz;  /* parking option */
     if (opaw != 1 && paw == 1)autopilot = 1-autopilot; /* car moves on its own */
     }

  if(opaw==2 && paw == 7){ /* restart, the most useful feature */
     deFault(); ohx=ohy=ohz=0; owx=owy=owz=0;
     }
        wnd=1;
	if(wnd){
      CAVEGetWandOrientation(azi,ele,rol);
      CAVEGetHead(hx,hy,hz);  /*position*/
      CAVEGetWand(wx,wy,wz);  /*position*/
      }
   else{  /* fix this for CAVE mouse flying */
      azi = .05*512; ele = .05*-512;
      }
   glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();

   glTranslatef(ohx-hx,ohy-hy,ohz-hz);
   if(mode==TURNMODE) glTranslatef(aff[12],aff[13],aff[14]);
   if(!friz){  /* if not parked */
      glRotatef(rol*(mode?torq:-torq),0.,0.,1.);
      glRotatef(ele*-torq, 1.,0.,0.);
      glRotatef(azi*-torq,0.,1.,0.);
      }
   if(joy&&((fabsf(CAVE_JOYSTICK_Y)>.02)||(fabsf(CAVE_JOYSTICK_X)>.02)))
       glTranslatef(CAVE_JOYSTICK_X*speed,0.,CAVE_JOYSTICK_Y*speed);

 if(opaw==2 && paw==2)glTranslatef(wx-owx, wy-owy, wz-owz); /*wand tractor */
 if(mode == TURNMODE) glTranslatef(-aff[12],-aff[13],-aff[14]);
   glMultMatrixf(aff);
   glGetFloatv(GL_MODELVIEW_MATRIX,aff);

  {int ii,jj; for(ii=0;ii<3;ii++){ luxx[ii]=0; /* calculite */
   for(jj=0;jj<3;jj++) luxx[ii] += aff[ii*4+jj]*lux[jj];} }

 if(mode==FLYMODE && !friz) /*make the star matrix rotate*/
      {
      glLoadIdentity();
      glRotatef(rol*(mode?torq:-torq),0.,0.,1.);
      glRotatef(ele*-torq,1.,0.,0.);
      glRotatef(azi*-torq,0.,1.,0.);
      glMultMatrixf(starmat);
      glGetFloatv(GL_MODELVIEW_MATRIX,starmat);
      }
   glPopMatrix();
   opaw=paw; ohx=hx; ohy=hy; ohz=hz;
             owx=wx; owy=wy; owz=wz;
   if(morph) autotymer(0);  /* advance autotymer */
/* audiofunc(); */
}

void drawcaveinit(void){   /* kludge? needed why ? */
   CAVEFar = 1000.;
   glClearColor(0.,0.,0.,0.);
   glClearDepth(1.);
}
#if 0
//This is how it works in the console
      float inv[16]; invert(inv, fram);  /* we want the inverse frame */
      glRotatef(180,0,1,0);  /*hack to correct the facing the wrong way problem*/
      glMultMatrixf(inv);
      glPushMatrix();
      glScalef(2,2,2);
      drawstars();
      glPopMatrix();
#endif

void drawcave(void){
  float hx,hy,hz;
  float inv[16], rinv[16];
  if(cockpit){  /* this could be neater */
    invert(inv, showchaseplane? fram2: fram); /* if there is a chaseplane, look from it */
    {int ii; FOR(ii,0,16)rinv[ii]=inv[ii]; }
    rinv[12]=rinv[13]=rinv[14]=0.;  /* just the rotation part of inv */
  }
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  drawstars(cockpit?rinv:starmat);   /* replace starmat by rot(fram^-1) */
  graffiti();
  CAVEGetHead(hx,hy,hz);
  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  glTranslatef(hx,hy,hz);  /* ancient and sacred translation to the head*/
  if(cockpit)glRotatef(180,0,1,0);  /* so we face the front of the coaster */
  glMultMatrixf(cockpit? inv : aff);
  glScalef(siz,siz,siz);
  glTranslatef(0.,0.,2.);  /* puts the head above the track */
  /* but it doesn't work right here. I don't like removing this */
  drawall();
  glPopMatrix();
}
#endif  /* CAVE */

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%end CAVE %%%%%%%%%%%%%%*/


//---------------------szg--------------------------

void normalize(float *vec){ int ii;
  float temp = vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2] + vec[3]*vec[3];
  if (temp == 1){ /* printf(" The vector is of unit length.\n"); */ }
  else{ temp = sqrt(temp); for(ii=0;ii<4;ii++)vec[ii] /= temp; }
}

void q2m(float *quat, float *mat){ /* quaternion yields a rotation */
  float w = quat[0], x = quat[1], y = quat[2], z = quat[3]; 
  float ww, xx, yy, zz, xy, xz, yz, wz, wy, wx;
  ww = w*w; xx = x*x; yy = y*y; zz = z*z; 
  xy = x*y; xz = x*z; yz = y*z; wz = w*z; wy = w*y; wx = w*x; xz = x*z;
  mat[0]=ww+xx-yy-zz;  mat[4]=2*xy - 2*wz;  mat[8]=2*wy + 2*xz;
  mat[1]=2*wz + 2*xy;  mat[5]=ww-xx+yy-zz;  mat[9]=2*yz - 2*wx;
  mat[2]=2*xz - 2*wy;  mat[6]=2*wx + 2*yz; mat[10]=ww-xx-yy+zz;
  mat[3] = mat[7] = mat[11] = mat[12] = mat[13] = mat[14] = 0;
  mat[15] = ww + xx + yy + zz; /* this is a check */
}

void m2q(float *mat, float *quat){  /* isometry yields a quaternion */ 
    float tmp; int ii, jj; 
    float tr[4] ={  /* pseudo traces */
      mat[0] + mat[5] + mat[10],
      mat[0] - mat[5] - mat[10],
     -mat[0] + mat[5] - mat[10],
     -mat[0] - mat[5] + mat[10]};
//  printf("trace  %f %f %f %f \n" , tr[0],tr[1],tr[2],tr[3]);  
    ii=1; jj=2;  /* find the largest tr */
    if(tr[0] > tr[1])   ii = 0; 
    if(tr[3] > tr[2])   jj = 3; 
    if(tr[jj] > tr[ii]) ii = jj;
//    printf("\n %i is the largest trace. \n", ii);
    tmp = 4*sqrt(tr[ii]+1)/2; 
    switch(ii){
    case 0: 
      quat[0] = tmp/4;
      quat[1] = (mat[6] - mat[9]) / tmp; 
      quat[2] = (mat[8] - mat[2]) / tmp; 
      quat[3] = (mat[1] - mat[4]) / tmp; 
      break; 
    case 1: 
      quat[0] = (mat[6] - mat[9]) / tmp;
      quat[1] = tmp/4; 
      quat[2] = (mat[1] + mat[4]) / tmp;
      quat[3] = (mat[8] + mat[2]) / tmp;
      break;
    case 2: 
      quat[0] = (mat[8] - mat[2]) / tmp;
      quat[1] = (mat[4] + mat[1]) / tmp; 
      quat[2] = tmp/4; 
      quat[3] = (mat[9] + mat[6]) / tmp; 
      break; 
    case 3:
      quat[0] = (mat[1] - mat[4]) / tmp;
      quat[1] = (mat[8] + mat[2]) / tmp; 
      quat[2] = (mat[9] + mat[6]) / tmp; 
      quat[3] = tmp/4;
      break; 
    }
}

void getmat(float *target, float *source){ 
     memcpy(target,source, 16*sizeof(float)); }

void multmatinv(float *T){  // of an Euclidean isometry 
   float N[16];  getmat(N,T);
   N[1] = T[4]; N[2] = T[8]; N[6] = T[9];
   N[4] = T[1]; N[8] = T[2]; N[9] = T[6];
   N[12]= -T[0]*T[12]-T[1]*T[13]-T[2]*T[14];
   N[13]= -T[4]*T[12]-T[5]*T[13]-T[6]*T[14];
   N[14]= -T[8]*T[12]-T[9]*T[13]-T[10]*T[14];
   glMultMatrixf(N);
}


void updateaffmat(float *oldhand, float *newhand){ 
  float dx=0, dy=0, dz=0 ; float dhand[16]; float quat[4];
  //we do arithmetic on an unused stack for convenience
  glMatrixMode(GL_TEXTURE); 
  glPushMatrix(); glLoadIdentity();
    multmatinv(newhand);
    glMultMatrixf(oldhand); 
    glGetFloatv(GL_TEXTURE_MATRIX, dhand);
  glPopMatrix(); //dhand= newhand^-1 * oldhand 

  dx=del*dhand[12]; dhand[12]=0;  //attenuation is about right
  dy=del*dhand[13]; dhand[13]=0; 
  dz=del*dhand[14]; dhand[14]=0; 
//use a bit of the translation in dhand and then strip translation 
  if(gesture){ //shrink dhand
    m2q(dhand,quat);
    quat[0]= 1-eps + eps*quat[0];
    quat[1] *=  eps; quat[2] *=  eps; quat[3] *=  eps;
    normalize(quat); 
    q2m(quat,dhand); 
  }

 if(mode==FLYMODE){
    glMatrixMode(GL_TEXTURE); glLoadIdentity();
    glMultMatrixf(dhand); //rotate the stars 
    glMultMatrixf(starmat);
    glGetFloatv(GL_TEXTURE_MATRIX,starmat);
 }

  glMatrixMode(GL_TEXTURE); glPushMatrix(); glLoadIdentity();
  if(mode==TURNMODE) glTranslatef(aff[12],aff[13],aff[14]);   
  glMultMatrixf(dhand);  //just the rotation of it
  glTranslatef(dx,dy,dz); 
  if(mode==TURNMODE) glTranslatef(-aff[12],-aff[13],-aff[14]);
  glMultMatrixf(aff); glGetFloatv(GL_TEXTURE_MATRIX,aff);
}

void readbuttons(arMasterSlaveFramework& fw){
 const int butn[6]= { fw.getButton(0), fw.getButton(1), fw.getButton(2), fw.getButton(3), fw.getButton(4), fw.getButton(5) };

 static int obutn[6]={0,0,0,0,0,0};
   
   if(obutn[4]==0 & butn[4] == 1)            ;  // {}-button
   if(obutn[1]==0 & butn[1] == 1) resethand=1;  //{B}-button
   if(obutn[2]==0 & butn[2] == 1) gesture=1-gesture; //{C}-button
   if(obutn[0]==0 & butn[0] == 1) changeScene=1;  //{A}-button
   if(obutn[3]==0 & butn[3] == 1) changeMol=1;  //{X}-button
   if(obutn[5]==0 & butn[5] == 1) deFault();    //{Z}-button
      
#ifdef DiVE
   if(obutn[1]==0 & butn[1] == 1) changeScene=1;  //{R}ed button
   if(obutn[3]==0 & butn[3] == 1) gesture=1-gesture; //{B}lue bbutton
   if(obutn[2]==0 & butn[2] == 1) resethand=1;  //{G}reen button
   if(obutn[4]==0 & butn[4] == 1) changeMol=1;  //{J}oy button
   if(obutn[5]==0 & butn[5] == 1) deFault();    //{T}rigger zap
   if(obutn[0]==0 & butn[0] == 1)               //{Y}ellow trail cycle
         {
         static int cnt=0; switch(++cnt){
         case 1 : trailmake=1; break;     //turn on trail
         case 2 : gap=6; break;           //erase surface
         case 3 : cnt= 0; trailmake=0; gap=gap0; break; //as before
         }}       
#endif
 memcpy(obutn,butn,6*sizeof(int));
} // end readbuttons


void gesturetrack(arMasterSlaveFramework& fw){
//the present logic is legacy stuff from zyspace, rationalize it!
  if (!fw.getMaster())return; //only the master should listen
//this is the gesture navigator
  static float  oldhand[16]; float handmat[16];

  arMatrix4 arhand=fw.getMatrix(1); getmat(handmat,arhand.v);
//geotrack(arhand.v) gets arhand.v and used as float* handmat
  static bool  reshand=true; if(resethand){reshand=true; resethand=false;}
//do we really need two flags?
  if(gesture){ updateaffmat(oldhand,handmat); }
  if(reshand){reshand=false; getmat(oldhand, handmat);}

//this is the Pape navigator
//this should really be done on the arithmetic (Texture) stack 
  glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();
  const float joy0 = fw.getAxis(0);  //joystick right-left 
  const float joy1 = fw.getAxis(1); //joystick fore-aft
  const float xx = arhand.v[8];
  const float yy = arhand.v[9];
  const float zz = arhand.v[10];
// printf(" xx %g yy %g zz %g \n",xx,yy,zz);

  const float da = .5*joy0; 
  const float ds = .1*joy1; 
//joystick translates in the wand direction
  glTranslatef(xx*ds, yy*ds, zz*ds);
  glRotatef(da, 0.,1.,0.);
  glMultMatrixf(aff);
  glGetFloatv(GL_MODELVIEW_MATRIX,aff);
  glPopMatrix();

  //this is an illogical place for this but maybe it works again

 {int ii,jj; for(ii=0;ii<3;ii++){ luxx[ii]=0; /* calculite */
              for(jj=0;jj<3;jj++) luxx[ii] += aff[ii*4+jj]*lux[jj];} }
					   
  readbuttons(fw);

} //end gesturetrack

bool cubeinit(arMasterSlaveFramework& fw, arSZGClient& ){

   ar_setNavMatrix(ar_translationMatrix(0,-5,0)); //initial position

//this is how Syzygy does shared memory.
   
  fw.addTransferField("aff", aff, AR_FLOAT, 16);
  fw.addTransferField("aff2", aff2, AR_FLOAT, 16);
  fw.addTransferField("starmat", starmat, AR_FLOAT, 16);

  fw.addTransferField("gesture", (&gesture), AR_INT, 1);
  fw.addTransferField("resethand", (&resethand), AR_INT, 1);

  fw.addTransferField("dfl2", (&dfl2), AR_INT, 1);
  fw.addTransferField("changeScene", (&changeScene), AR_INT, 1);
  fw.addTransferField("maxMol", (&maxMol), AR_INT, 1);
  fw.addTransferField("changeMol", (&changeMol), AR_INT, 1);
  fw.addTransferField("sceneType", (&sceneType), AR_INT, 1);
  fw.addTransferField("openFileNum", (&openFileNum), AR_INT, 1);
  fw.addTransferField("movieCounter", (&movieCounter), AR_INT, 1);
  fw.addTransferField("selectChains", (&selectChains), AR_INT, 1);
  fw.addTransferField("whichChain", (&whichChain), AR_INT, 1);
  fw.addTransferField("showHydrophob", (&showHydrophob), AR_INT, 1);
  fw.addTransferField("backBone", (&backBone), AR_INT, 1);
 
  fw.addTransferField("sizz", (&sizz), AR_FLOAT, 1);
  fw.addTransferField("mode", (&mode), AR_INT, 1);
  fw.addTransferField("luxx", luxx, AR_FLOAT, 3);
  fw.addTransferField("mauspaw", &mauspaw, AR_INT, 1);
  fw.addTransferField("th0", (&th0), AR_INT, 1);
  fw.addTransferField("th1", (&th1), AR_INT, 1);
  fw.addTransferField("ta0", (&ta0), AR_INT, 1);
  fw.addTransferField("ta1", (&ta1), AR_INT, 1);
  fw.addTransferField("one", (&one), AR_FLOAT, 1);




  fw.addTransferField("gap", (&gap), AR_FLOAT, 1);
  fw.addTransferField("morph", (&morph), AR_INT, 1);
//#endif
  return true;
}

void clefs(arMasterSlaveFramework& , unsigned char key, int x, int y){ 
  keyboard(key,x,y);
}


void postexchange(arMasterSlaveFramework& fw){ /* console idle */ 
 // cerr << "postexchange.\n";
	if(dfl2){dfl2=0;deFault();reSet();}
  // cerr << "postexchange2.\n";
	if(changeScene){
		changeScene=0;
		if(sceneType<4){sceneType=sceneType+1;} 
		else {sceneType=0;};
		reSet();
	}
// cerr << "sceneType: " << sceneType << endl;
 //  cerr << "postexchange3.\n";
	if(changeMol){
		changeMol=0;
		if(openFileNum<(maxMol-1)){openFileNum=openFileNum+1;} 
		else {openFileNum=0;};
		openpdbfile();
	}
  // cerr << "postexchange4.\n";
}


void cubedraw(arMasterSlaveFramework& fw){
  // cerr << "draw.\n";
  glPushMatrix();
  glTranslatef(0,5,0);
  glScalef(sizz,sizz,sizz);
  // cerr << "draw2.\n";
  // drawstars(starmat);
  // cerr << "draw3.\n";
  glPopMatrix();
  glMultMatrixf(aff);
  glScalef(sizz,sizz,sizz);
  // dawn(1);  doesnothing here
  drawall();
  dawn(1);
   //  cerr << "draw4.\n";
  if(mov) movie();
  // cerr << "done draw.\n";
  }

#if 0
void postit(void){     
  char buf[256]; /* console postit are done differently from cave */
#define  LABEL2(x,y,W,u) {sprintf(buf,(W),(u));char2wall(x,y,0.,buf);}
  glMatrixMode(GL_PROJECTION); glPushMatrix(); glLoadIdentity();
  gluOrtho2D(0,3000,0,3000);
  glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity();
      /*bull's eye*/
      if(mode==TURNMODE) glColor3f(0x22/255.,0x88/255.,0xdd/255.);
        else      glColor3f(.8,.8,.8);
      LABEL2(1500,1500,"%s","o");
      /* writings */
      if(mode==TURNMODE) glColor3f(1.,0.,1.);
        else      glColor3f(1.,1.,0.);
      LABEL2(80,80,"%4.1f fps",speedometer());

      if(eyes==cube){
#if 0
      LABEL2(80,2140,"moving %d",moving);
      LABEL2(80,2000,"(t)railmake %d",trailmake);
      LABEL2(80,1860,"trailtie %d",trailtie);
#endif
  
      LABEL2(80,2840,\
     "CUBE  {0} {1}recalibrate {2}gesture {3}T-key {4}I-key {5}zap )",
			           mode?"fly":"turn");
#ifdef DiVE
     "DiVE  {R}morph {Y}crumbs {G}recalibrate {B}gesture %s{J}mode {T}zap )",
			           mode?"fly":"turn");
#endif
      LABEL2(10,10,"%s","illiSnail 21JAN06 by Francis, Bourd, Hartman, Chappell (C) 1995..2005, U Illinois.");
      }else{//if console
      LABEL2(80,2840,\
      "(ESC)ape (V)Binoc (MAUS2)Fore ) %s(M)ode (H)omotopy (W)riting", mode?"TURN":"FLY");
      LABEL2(10,10,"illiSnail 06nov05  \
   by Francis, Bourd, Hartman & Chappell, U Illinois, 1995..1997 %s","");
      LABEL2(80,2770,"(N)ose   %0.3f",nose);
      LABEL2(80,2700,"(S)peed  %0.4f",speed);
      LABEL2(80,2630," tor(Q) %0.4f",torq);
      LABEL2(80,2560,"near clipper %g", mysiz*focal);
      LABEL2(80,2490,"f(O)cal factor %g",focal);
      LABEL2(80,2420,"my s(I)ze %.2g",mysiz);
      LABEL2(80,2280,"(Z)ap %s","");
      LABEL2(80,2210,"(G)ap %.2g",gap);
      } //endif console
    glPopMatrix();
    glMatrixMode(GL_PROJECTION); glPopMatrix();
} 
#endif

void overlay(arMasterSlaveFramework&){ messages2(); }


//----------------------szg------------------------------------


int main(int argc, char **argv){  
   eyes=cube;
//   arguments(argc,argv); 
   if(eyes==cave)
#ifdef CAVE
       CAVEConfigure(&argc, argv, NULL)
#endif
       ;
   getmem();
   dataprep();
 
   if(eyes==cave){
#ifdef CAVE
       CAVEInit(); /* Each wall is (part of) a forked process from here on.*/ 
       CAVEFrameFunction(cavetrack,0);  
       /*is restricted to MasterWall, so once per frame */
       CAVEInitApplication(drawcaveinit, 0); 
       CAVEDisplay(drawcave, 0);
       while(!CAVEgetbutton(CAVE_ESCKEY))
         {
         cavekeybo();  /* is asynchronous from display processes */ 
	 }
       audioclean();
       CAVEExit();
#endif
       ;}
     else  if(eyes==console){ /*console main*/
       glutInit(&argc, argv);
       glutInitDisplayMode(GLUT_DOUBLE|GLUT_DEPTH);
       switch(win)
         { case 0: break;
           case 1: glutInitWindowSize(640, 480);
                   glutInitWindowPosition(0,1024-480);
                   break;
           case 2: glutInitWindowPosition(0,0);
	           break;
         }
       glutCreateWindow("<* illiSkel in C/OpenGL/GLUT *>");
       if(win==2) glutFullScreen();
       glEnable(GL_DEPTH_TEST);
       glutDisplayFunc(drawcons);       
       glutKeyboardFunc(keyboard);
       glutSpecialFunc(special_keybo);
       glutMouseFunc(mousepushed);
       glutMotionFunc(mousemoved);       
       glutPassiveMotionFunc(mousemoved); 
       glutReshapeFunc(reshaped);
       glutIdleFunc(idle);             
       glutMainLoop();
       }
     else if(eyes==cube){ //syzygy main
       //fff = new arMasterSlaveFramework;  which is new?
       arMasterSlaveFramework* fff = new arMasterSlaveFramework();
       
       if(!fff->init(argc,argv)) return 1;
       fff->setStartCallback(cubeinit); 
   //     fff->setInitCallback(cubeinit); 
   //    fff->setPreExchangeCallback(cubetrack);
       fff->setPreExchangeCallback(gesturetrack);
       fff->setDrawCallback(cubedraw);
       fff->setPostExchangeCallback(postexchange);
       fff->setKeyboardCallback(clefs);
       fff->setOverlayCallback(overlay);
         // framework.setClipPlanes( .1, 1000 );
         // framework.setEyeSpacing(6/(12*2.54));
         //return framework.start() ? 0 : 1;
       return fff->start() ? 0:1;
     }
}