/*
 *   povchem.c; July 12, 1996.
 */
 
char *VERSION = "1.00";

/*
 *   Copyright (C) 1995-1996 by Paul Thiessen (paul@ludwig.scs.uiuc.edu)
 *
 *   This program is shareware and may be used freely for non-profit work.
 *   It is not free for commercial use, however; see the legal info in the
 *   Manual page (below).
 *
 *   By using this program you are agreeing to abide by all applicable
 *   copyright laws. You also are agreeing that the author cannot be held
 *   liable for any losses incurred by the use of this program.
 *
 *   For more on legal rights and authorship, and for instructions on how 
 *   to use the program, see
 *
 *            http://ludwig.scs.uiuc.edu/~paul/Manual.html
 */

/*
 * To do someday, maybe:
 *   - make clever color assignment manipulator
 *   - make dynamic PDB tag constructor
 *   - put in bond calculator?
 */

/*------------ THESE ARE THE ONLY USER-CONFIGURABLE OPTIONS -------------*/

/*
 * The maximum number of elements in the Hydrogen bond donors and acceptors
 * lists.
 */
#define MAX_H_BOND_ELEMENTS 10

/*
 * Basically these just define the floating point type (F_TYPE), the number
 * of significant figures (roughly) to which that type is accurate, and
 * the names of various mathematical functions. All of these probably won't
 * need to be changed, but if memory is extremely tight and you're processing
 * really big molecules, you can set F_TYPE to 'float' instead of 'double'
 * and on most machines this will save some space.
 */
#define F_TYPE double
#define SF 15
#define SQRT sqrt
#define COS cos
#define SIN sin
#define ACOS acos
#define ASIN asin
#define POW pow


/*--------------- DO NOT MODIFY ANYTHING BELOW THIS LINE! ----------------*/

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

/* Basic models */
enum { CPK=1, BALLnSTICK, CYLINDER };

/* Coloring schemes */
enum { BY_ELEMENT=1, BY_PDB, ALL_THE_SAME, BY_PDB_SEPARATE };

/* BY_PDB options */
enum { PDB_RES=1, PDB_CH, PDB_RES_CH_SEQ, PDB_RES_ELEM };

/* How sticks/cylinders change color */
enum { SPLIT_EQUAL=1, SPLIT_PROPORTIONAL, SPLIT_SHADED };

/* Size of balls in ball and stick model */
enum { CONSTANT=1, PROPORTIONAL };

/* Image orientation */
enum { LANDSCAPE=1, PORTRAIT };

/* Background colors */
enum { BLACK=1, BLUE, WHITE };

/* Views */
enum { BEST_GUESS=1, STANDARD, USER_DEFINED };

#define FALSE 0
#define TRUE 1

#define U_INT unsigned int

/*
 * Configuration file options. Default variables are in all caps, next
 * to their respective "real" variables.
 */

char *PERIODIC_TABLE_FILE=      NULL;

char *POVRAY=           NULL;         /* text variable assigned later - */
U_INT POV_WIDTH=        0;            /* after parsing config file to   */
char *POV_OPTIONS=      NULL;         /* avoid string memory waste.     */
U_INT RUN_POVRAY=       FALSE;

char *VIEWER=           NULL;

char *FINISH=                   NULL;
char *FINISH_DECLARATION=       NULL;

U_INT  DEFAULT_MODEL=                   CYLINDER,       model;
U_INT  DEFAULT_COLORMODEL=              BY_ELEMENT,     cylinderColorModel,
                                                        sphereColorModel;
U_INT  DEFAULT_PDB_FIELDS=              PDB_RES_CH_SEQ, sphereTagType=0;
F_TYPE DEFAULT_CYLINDER_RADIUS=         0.25,           cylinderRadius;
U_INT  DEFAULT_RADIUS_TYPE=             PROPORTIONAL,   ballSizeModel;
F_TYPE DEFAULT_PROPORTIONALITY=         0.20;
F_TYPE DEFAULT_BALL_RADIUS=             0.3,            ballRadius;
U_INT  DEFAULT_THICKNESS=               CONSTANT,       cylinderSizeModel;
F_TYPE DEFAULT_STICK_RADIUS=            0.1;
F_TYPE DOUBLE_FACTOR=                   1.5;
F_TYPE TRIPLE_FACTOR=                   2.0;
F_TYPE HIGH_ORDER_FACTOR=               2.5;
U_INT  DEFAULT_STICKS_LIKE_BALLS=       FALSE;
U_INT  DEFAULT_STICK_COLORMODEL=        BY_PDB;
U_INT  DEFAULT_STICK_PDB_FIELDS=        PDB_RES,        stickTagType=0;
U_INT  DEFAULT_BOND_SPLIT=              SPLIT_PROPORTIONAL, stickSplit;
F_TYPE POWER=                           2.0;
F_TYPE DEFAULT_BALL_RED=                0.6;
F_TYPE DEFAULT_BALL_GREEN=              0.2;
F_TYPE DEFAULT_BALL_BLUE=               0.9;
F_TYPE DEFAULT_STICK_RED=               0.2;
F_TYPE DEFAULT_STICK_GREEN=             0.8;
F_TYPE DEFAULT_STICK_BLUE=              0.5;
U_INT  DEFAULT_FIND_H_BONDS=            FALSE;
U_INT  DEFAULT_H_BOND_DASHES=           3,              HBondDashes;
F_TYPE DEFAULT_H_BOND_RADIUS=           0.1,            HBondRadius;
char  *H_BOND_DONORS=                   NULL;
char  *H_BOND_ACCEPTORS=                NULL;
F_TYPE DEFAULT_MINDIST=                 1.0,            HBondMinDist;
F_TYPE DEFAULT_MAXDIST=                 3.1,            HBondMaxDist;
F_TYPE DEFAULT_MIN_DHA_ANGLE=           140.0,          HBondMinDHA;
F_TYPE DEFAULT_MIN_HAR_ANGLE=           60.0,           HBondMinHAR;
F_TYPE DEFAULT_H_BOND_RED=              1.0;
F_TYPE DEFAULT_H_BOND_GREEN=            1.0;
F_TYPE DEFAULT_H_BOND_BLUE=             0.0;
F_TYPE DEFAULT_BACKGROUND_RED=          0.0;
F_TYPE DEFAULT_BACKGROUND_GREEN=        0.0;
F_TYPE DEFAULT_BACKGROUND_BLUE=         0.0;
U_INT  DEFAULT_ORIENTATION=             LANDSCAPE,      imageOrientation;
U_INT  DEFAULT_VIEW=                    STANDARD,       viewType;
U_INT  DEFAULT_MAKE_STEREO_PAIR=        FALSE;
F_TYPE STEREO_SEPARATION_FACTOR=        0.10;

F_TYPE CAMERA_X=        0.0;
F_TYPE CAMERA_Y=        0.0;
F_TYPE CAMERA_Z=        -20.0;
F_TYPE LOOK_AT_X=       0.0;
F_TYPE LOOK_AT_Y=       0.0;
F_TYPE LOOK_AT_Z=       0.0;
F_TYPE SKY_X=           0.0;
F_TYPE SKY_Y=           1.0;
F_TYPE SKY_Z=           0.0;
F_TYPE DIRECTION=       2.0;
F_TYPE WINDOW_UP=       1.0;
F_TYPE WINDOW_RIGHT=    1.333;

U_INT MAX_BONDS_PER_ATOM=       5;

/*----------- end configuration file options ---------------*/


/*---------- macros -----------*/

#define VDW(INDEX) ((PeriodicTable[(Atoms[INDEX].Z)]).vdW)
#define NAME(INDEX) ((PeriodicTable[(Atoms[INDEX].Z)]).name)

#define ZNUM(INDEX) (Atoms[INDEX].Z)
#define SERNUM(INDEX) (Atoms[INDEX].serNum)
#define TAG(INDEX) (Atoms[INDEX].tag)
#define X(INDEX) (Atoms[INDEX].x)
#define Y(INDEX) (Atoms[INDEX].y)
#define Z(INDEX) (Atoms[INDEX].z)

#define RED(Z) ((PeriodicTable[Z]).r)
#define GREEN(Z) ((PeriodicTable[Z]).g)
#define BLUE(Z) ((PeriodicTable[Z]).b)

#define DEGREES(rad) ((180.0/M_PI)*rad)
#define RADIANS(deg) ((M_PI/180.0)*deg)

unsigned long memAlloced=0;

#define STRNCPY0(buf,pos,n) strncpy(buf,pos,n); buf[n]='\0'

#define BOMB fclose(pdb); fclose(pov); fclose(inc); exit(EXIT_FAILURE)


/*------------------- global variables and type def's ---------------*/

typedef struct { F_TYPE x, y, z;
                 U_INT serNum;
                 unsigned char Z;
                 char *tag; } Atom;
                 
Atom *Atoms;
U_INT nAtoms;

U_INT *BondsTo;
unsigned char *BondOrders;
#define BoundTo(bond,index) (BondsTo[(bond)*MAX_BONDS_PER_ATOM+(index)])
#define BondOrder(bond,index) (BondOrders[(bond)*MAX_BONDS_PER_ATOM+(index)])

U_INT HBondDonors[MAX_H_BOND_ELEMENTS], nDonors;
U_INT HBondAcceptors[MAX_H_BOND_ELEMENTS], nAcceptors;

enum { ZERO=0, SINGLE, DOUBLE, TRIPLE, HIGH_ORDER, H_BOND }; /* bond orders */

U_INT nOrders[6];
char *orderTypes[]={"Zero_","Single_","Double_",
                    "Triple_","High_Order_","H_"};

typedef struct { char *name;
                 char symbol[3];
                 F_TYPE vdW, r, g, b;
                 unsigned char isPresent; } Element;

#define NUM_ELEMENTS 111
Element PeriodicTable[NUM_ELEMENTS+1];
                
FILE *pdb, *inc, *pov;

U_INT isInteractive=TRUE;

F_TYPE aveX=0.0, aveY=0.0, aveZ=0.0;
F_TYPE minX=0.0, maxX=0.0, minY=0.0, maxY=0.0, minZ=0.0, maxZ=0.0;

typedef struct { F_TYPE x, y, z; } Vector;
Vector unitRight, Light, LeftEye, RightEye;

typedef struct { F_TYPE r, g, b; } Color;
Color Background, allSphereColor, allCylinderColor, HBondColor;

typedef struct { char *name;
                 F_TYPE r, g, b;
                 void *next; } tagType;

tagType *firstSphereTag=NULL, *firstCylinderTag=NULL;
U_INT nSphereTags=0, nCylinderTags=0;

typedef struct { U_INT serial;
                 U_INT Z;
                 char resName[4];
                 char chainID;
                 U_INT resSeq;
                 F_TYPE x, y, z; } PDBatomInfoType;

#define FIND_SERIAL  0x01
#define FIND_Z       0x02
#define FIND_RESNAME 0x04
#define FIND_CHAINID 0x08
#define FIND_RESSEQ  0x10
#define FIND_COORDS  0x20
#define FIND_ALL (FIND_SERIAL | FIND_Z | FIND_RESNAME | FIND_CHAINID | \
                  FIND_RESSEQ | FIND_COORDS)

char **cylinderEndTags;

enum { MONO=1, LEFT_EYE, RIGHT_EYE };
char stereoPairCreated=FALSE;

char *incName, *povName, *leftName=NULL, *rightName=NULL;


/*-------- return a U_INT rounded from an F_TYPE ----------*/

U_INT Round(F_TYPE fVal)
{
  double iVal, frac;
  
  frac=modf((double)fVal,&iVal);
  if (frac>=0.5) iVal+=1.0;
  return ((U_INT)iVal);
}
  

/*---- returns TRUE if a string contains only ' ', '\t', and '\n' ------*/

int isBlank(char *line)
{
  char *pos;
  int blank=TRUE;

  if (line) {
    for (pos=line; *pos!='\0' && blank==TRUE; pos++) {
      if (*pos!=' ' && *pos!='\t' && *pos!='\n') blank=FALSE;
    }
  } else
    blank=FALSE;
  return blank;
}

/*--- converts a string to lower case, with first letter upper case ---*/

char *str_tolower(char *line)
{
  char *pos;

  if (line && *line) {
      *line=toupper(*line);
      for (pos=line+1; *pos!='\0'; pos++)
        *pos=tolower(*pos);
  }
  return line;
}

/*---- makes a new string, keeping track of memory allocation ----*/

char *New_String(char **new, const char *old)
{
  if ( ((*new)=strdup(old))==NULL ) {
    puts("FATAL: out of memory during strdup.");
    exit(EXIT_FAILURE);
  } else
    memAlloced+=strlen(*new)+1;
  return *new;
}


/*------ Parsing functions - for periodic table and config file -------*/

/*------- parses the hydrogen bond donor/acceptor lists --------*/

U_INT Determine_Z_From_Symbol(const char *symbol)
{
  U_INT Z;

  for (Z=1; Z<NUM_ELEMENTS; Z++) {
    if (strcmp(PeriodicTable[Z].symbol,symbol)==0) return Z;
  }
  return 0;
}

void Make_H_Bonder_List(const char *list)
{
  U_INT Z, ZNums[MAX_H_BOND_ELEMENTS], count=0;
  char copy[128], *symbol;

  strcpy(copy,list);
  for (symbol=str_tolower(strtok(copy," \t,"));
       symbol!=NULL;
       symbol=str_tolower(strtok(NULL," \t,"))  ) {
    if (count==MAX_H_BOND_ELEMENTS) {
      printf("Warning: too many elements in hydrogen bond %s list.\n",
             (list==H_BOND_DONORS) ? "donor" : "acceptor");
      printf(
        "Elements beyond MAX_H_BOND_ELEMENTS, currently %u, will be ignored.\n",
             MAX_H_BOND_ELEMENTS);
      break;
    }
    if ( (ZNums[count]=Determine_Z_From_Symbol(symbol))==0 ) {
      printf("Element '%s' in hydrogen bond %s list not recognized,\n",
        symbol, (list==H_BOND_DONORS) ? "donor" : "acceptor");
      puts("so it will be ignored.\n");
    } else
      count++;
  }
  if (list==H_BOND_DONORS) {
    nDonors=count;
    for (Z=0;Z<count;Z++) HBondDonors[Z]=ZNums[Z];
  } else {
    nAcceptors=count;
    for (Z=0;Z<count;Z++) HBondAcceptors[Z]=ZNums[Z];
  }
/*
  printf("Found h-bonder list: ");
  for (Z=0;Z<count;Z++) printf("%s ",PeriodicTable[ZNums[Z]].symbol);
  puts("");
*/
  return;
}


/*
 * This function assigns default values to many variables, and, if present,
 * the configuration file is read to assign new default values.
 */

typedef struct { const char *name;
                 void *variable; } configVarMap;

#define ERROR { error=TRUE; goto parsing_error; }

void Assign_Defaults(char **fileName)
{
  FILE *config;
  char line[256], *option, *value, *tmp;
  U_INT lineNum=0, error, i;
  int iVal;
  void *var;
  F_TYPE fVal;

/*
 * This is a mapping of all the option variables in the config file, to their
 * respective internal program variables.
 */

  configVarMap configVars[]= {
    { "Periodictable",          &PERIODIC_TABLE_FILE            },
    { "Povray",                 &POVRAY                         },
    { "Povwidth",               &POV_WIDTH                      },
    { "Povoptions",             &POV_OPTIONS                    },
    { "Runpovray",              &RUN_POVRAY                     },
    { "Finish",                 &FINISH                         },
    { "Userfinish",             &FINISH_DECLARATION             },
    { "Viewer",                 &VIEWER                         },
    { "Model",                  &DEFAULT_MODEL                  },
    { "Coloring",               &DEFAULT_COLORMODEL             },
    { "Pdbfields",              &DEFAULT_PDB_FIELDS             },
    { "Cylinderradius",         &DEFAULT_CYLINDER_RADIUS        },
    { "Ballradiustype",         &DEFAULT_RADIUS_TYPE            },
    { "Proportion",             &DEFAULT_PROPORTIONALITY        },
    { "Ballradius",             &DEFAULT_BALL_RADIUS            },
    { "Orderthickness",         &DEFAULT_THICKNESS              },
    { "Stickradius",            &DEFAULT_STICK_RADIUS           },
    { "Doublefactor",           &DOUBLE_FACTOR                  },
    { "Triplefactor",           &TRIPLE_FACTOR                  },
    { "Highorderfactor",        &HIGH_ORDER_FACTOR              },
    { "Colorstickslikeballs",   &DEFAULT_STICKS_LIKE_BALLS      },
    { "Stickcoloring",          &DEFAULT_STICK_COLORMODEL       },
    { "Stickpdbfields",         &DEFAULT_STICK_PDB_FIELDS       },
    { "Colorsplit",             &DEFAULT_BOND_SPLIT             },
    { "Power",                  &POWER                          },
    { "Findhbonds",             &DEFAULT_FIND_H_BONDS           },
    { "Hbonddashes",            &DEFAULT_H_BOND_DASHES          },
    { "Hbonddonors",            &H_BOND_DONORS                  },
    { "Hbondacceptors",         &H_BOND_ACCEPTORS               },
    { "Hbondradius",            &DEFAULT_H_BOND_RADIUS          },
    { "Mindist",                &DEFAULT_MINDIST                },
    { "Maxdist",                &DEFAULT_MAXDIST                },
    { "Mindhaangle",            &DEFAULT_MIN_DHA_ANGLE          },
    { "Minharangle",            &DEFAULT_MIN_HAR_ANGLE          },
    { "Hbondred",               &DEFAULT_H_BOND_RED             },
    { "Hbondgreen",             &DEFAULT_H_BOND_GREEN           },
    { "Hbondblue",              &DEFAULT_H_BOND_BLUE            },
    { "Allatomred",             &DEFAULT_BALL_RED               },
    { "Allatomgreen",           &DEFAULT_BALL_GREEN             },
    { "Allatomblue",            &DEFAULT_BALL_BLUE              },
    { "Allbondred",             &DEFAULT_STICK_RED              },
    { "Allbondgreen",           &DEFAULT_STICK_GREEN            },
    { "Allbondblue",            &DEFAULT_STICK_BLUE             },
    { "Backgroundred",          &DEFAULT_BACKGROUND_RED         },
    { "Backgroundgreen",        &DEFAULT_BACKGROUND_GREEN       },
    { "Backgroundblue",         &DEFAULT_BACKGROUND_BLUE        },
    { "Orientation",            &DEFAULT_ORIENTATION            },
    { "View",                   &DEFAULT_VIEW                   },
    { "Makestereopair",         &DEFAULT_MAKE_STEREO_PAIR       },
    { "Stereoseparationfactor", &STEREO_SEPARATION_FACTOR       },
    { "Camerax",                &CAMERA_X                       },
    { "Cameray",                &CAMERA_Y                       },
    { "Cameraz",                &CAMERA_Z                       },
    { "Lookatx",                &LOOK_AT_X                      },
    { "Lookaty",                &LOOK_AT_Y                      },
    { "Lookatz",                &LOOK_AT_Z                      },
    { "Skyx",                   &SKY_X                          },
    { "Skyy",                   &SKY_Y                          },
    { "Skyz",                   &SKY_Z                          },
    { "Windowdistance",         &DIRECTION                      },
    { "Windowheight",           &WINDOW_UP                      },
    { "Windowwidth",            &WINDOW_RIGHT                   },
    { "Maxbondsperatom",        &MAX_BONDS_PER_ATOM             },
    { "end-of-list",            NULL                            }
  };

/*
 * The "default defaults" in the source code can be superceded by the
 * defaults in the config file.
 *
 * Options are parsed in the order present, case insensitive, and not all
 * variables need be assigned, since they have "default defaults" anyway.
 *
 */
  
  if (*fileName==NULL) {
    New_String(fileName,"povchem.cfg");
    if ((config=fopen(*fileName,"r"))==NULL) {
      puts("Can't read the configuration file \"povchem.cfg\"!");
      perror("Reason");
      puts("Assuming all defaults.\n");
      goto no_config;
    }
  } else {
    if ((config=fopen(*fileName,"r"))==NULL) {
      printf("Can't read the configuration file \"%s\"!\n",*fileName);
      perror("Reason");
      exit(EXIT_FAILURE);
    }
  }

  while (fgets(line,256,config)!=NULL) {
    lineNum++;
    error=FALSE;

    if (line[0]!='#' && !isBlank(line)) {
      if (!strchr(line,'=')) ERROR;
      option=str_tolower(strtok(line," \t="));
      if (!option) ERROR;

/*----------- determine which option the current line contains ----------*/

      var=NULL;
      for ( i=0; configVars[i].variable ; i++ ) {
        if (strcmp(option,configVars[i].name)==0) {
          var=configVars[i].variable;
          break;
        }
      }
      if (!var) ERROR;

/*
 * The value is everything from the first non-space after the '=', up to
 * the last non-space. If the line is blank after the '=' then that option
 * is ignored and assumes its default value.
 */

      if (!(value=strtok(NULL,"\n"))) ERROR;
      for (;
           *value!='\0' && (*value==' ' || *value=='=' || *value=='\t') ;
           value++ ) {}
      if (*value=='\0') goto parsing_error;
      for (tmp=value+strlen(value)-1;
           *tmp==' ' || *tmp=='\t';
           tmp-- ) { *tmp='\0'; }


      if (var==&PERIODIC_TABLE_FILE || var==&POVRAY || var==&FINISH || 
          var==&VIEWER || var==&H_BOND_DONORS || var==&H_BOND_ACCEPTORS) {
        New_String(var,value);
      }

      else if (var==&POV_OPTIONS) {
        if (POVRAY) New_String(var,value);
      }

      else if (var==&POV_WIDTH) {
        if (POVRAY) {
          if ( (POV_WIDTH=atoi(value))==0 ) {
            POVRAY=NULL;
            ERROR;
          }
        }
      }

      else if (var==&RUN_POVRAY) {
        if (POVRAY) {
          if (strcmp(str_tolower(value),"Yes")==0)
            RUN_POVRAY=TRUE;
          else if (strcmp(value,"No")==0)
            RUN_POVRAY=FALSE;
          else
            ERROR;
        } else
          RUN_POVRAY=FALSE;
      }

      else if (var==&FINISH_DECLARATION) {
        if ((FINISH_DECLARATION=(char *)malloc(strlen(value)+1))==NULL) {
          fclose(config);
          puts("\nFatal: Not enough memory to allocate user finish!\n\n");
          exit(EXIT_FAILURE);
        }
        sprintf(FINISH_DECLARATION,"%s\n",value);
        do {
          if (fgets(line,256,config)==NULL || 
              (tmp=(char *)
                  malloc(strlen(FINISH_DECLARATION)+strlen(line)+1) )==NULL) {
            puts("Fatal: UserFinish has no terminating '}' line,");
            puts("or is too long to fit in available memory!\n\n");
            exit(EXIT_FAILURE);
          } else {
            strcpy(tmp,FINISH_DECLARATION);
            free(FINISH_DECLARATION);
            FINISH_DECLARATION=tmp;
          }
          lineNum++;
          strcat(FINISH_DECLARATION,line);
        } while (*line!='}');
        strcat(FINISH_DECLARATION,"\n");
      }
      
      else if (var==&DEFAULT_MODEL) {
        if (strcmp(str_tolower(value),"Cpk")==0)
          DEFAULT_MODEL=CPK;
        else if (strcmp(value,"Ballandstick")==0)
          DEFAULT_MODEL=BALLnSTICK;
        else if (strcmp(value,"Cylinder")==0)
          DEFAULT_MODEL=CYLINDER;
        else
          ERROR;
      }

      else if (var==&DEFAULT_COLORMODEL || var==&DEFAULT_STICK_COLORMODEL) {
        if (strcmp(str_tolower(value),"Byelement")==0)
          *((U_INT *)var)=BY_ELEMENT;
        else if (strcmp(value,"Bypdb")==0)
          *((U_INT *)var)=BY_PDB;
        else if (strcmp(value,"Allthesame")==0)
          *((U_INT *)var)=ALL_THE_SAME;
        else
          ERROR;
      }
      
      else if (var==&DEFAULT_PDB_FIELDS || var==&DEFAULT_STICK_PDB_FIELDS) {
        if (strcmp(str_tolower(value),"Residue")==0)
          *((U_INT *)var)=PDB_RES;
        else if (strcmp(value,"Chain")==0)
          *((U_INT *)var)=PDB_CH;
        else if (strcmp(value,"Reschseq")==0)
          *((U_INT *)var)=PDB_RES_CH_SEQ;
        else if (strcmp(value,"Reselem")==0)
          *((U_INT *)var)=PDB_RES_ELEM;
        else
          ERROR;
      }
      
      else if (var==&DEFAULT_CYLINDER_RADIUS || var==&DEFAULT_PROPORTIONALITY ||
               var==&DEFAULT_BALL_RADIUS || var==&DEFAULT_STICK_RADIUS ||
               var==&DOUBLE_FACTOR || var==&TRIPLE_FACTOR || var==&POWER ||
               var==&HIGH_ORDER_FACTOR || var==&DEFAULT_H_BOND_RADIUS ||
               var==&DEFAULT_MINDIST || var==&DEFAULT_MAXDIST ||
               var==&DEFAULT_MIN_DHA_ANGLE || var==&DEFAULT_MIN_HAR_ANGLE ||
               var==&STEREO_SEPARATION_FACTOR || 
               var==&DIRECTION || var==&WINDOW_UP || var==&WINDOW_RIGHT) {
        if ((fVal=atof(value))<=0.0) ERROR;
        *((F_TYPE *)var)=fVal;
      }
      
      else if (var==&DEFAULT_RADIUS_TYPE || var==&DEFAULT_THICKNESS) {
        if (strcmp(str_tolower(value),"Proportional")==0)
          *((U_INT *)var)=PROPORTIONAL;
        else if (strcmp(value,"Constant")==0)
          *((U_INT *)var)=CONSTANT;
        else
          ERROR;
      }
      
      else if (var==&DEFAULT_STICKS_LIKE_BALLS || var==&DEFAULT_FIND_H_BONDS ||
               var==&DEFAULT_MAKE_STEREO_PAIR) {
        if (strcmp(str_tolower(value),"Yes")==0)
          *((U_INT *)var)=TRUE;
        else if (strcmp(value,"No")==0)
          *((U_INT *)var)=FALSE;
        else
          ERROR;
      }

      else if (var==&DEFAULT_BOND_SPLIT) {
        if (strcmp(str_tolower(value),"Equal")==0)
          DEFAULT_BOND_SPLIT=SPLIT_EQUAL;
        else if (strcmp(value,"Proportional")==0)
          DEFAULT_BOND_SPLIT=SPLIT_PROPORTIONAL;
        else if (strcmp(value,"Shaded")==0)
          DEFAULT_BOND_SPLIT=SPLIT_SHADED;
        else
          ERROR;
      }

      else if (var==&DEFAULT_H_BOND_DASHES || var==&MAX_BONDS_PER_ATOM) {
        if ( (iVal=atoi(value)) <=0) ERROR;
        *((U_INT *)var)=iVal;
      }

      else if (var==&DEFAULT_H_BOND_RED || var==&DEFAULT_H_BOND_GREEN ||
               var==&DEFAULT_H_BOND_BLUE || var==&DEFAULT_BACKGROUND_GREEN ||
               var==&DEFAULT_BACKGROUND_RED || var==&DEFAULT_BACKGROUND_BLUE ||
               var==&DEFAULT_BALL_RED || var==&DEFAULT_BALL_GREEN ||
               var==&DEFAULT_BALL_BLUE || var==&DEFAULT_STICK_RED ||
               var==&DEFAULT_STICK_GREEN || var==&DEFAULT_STICK_BLUE) {
        if ( (fVal=atof(value))<0.0 || fVal>1.0) ERROR;
        *((F_TYPE *)var)=fVal;
      }

      else if (var==&DEFAULT_ORIENTATION) {
        if (strcmp(str_tolower(value),"Portrait")==0)
          DEFAULT_ORIENTATION=PORTRAIT;
        else if (strcmp(value,"Landscape")==0)
          DEFAULT_ORIENTATION=LANDSCAPE;
        else
          ERROR;
      }

      else if (var==&DEFAULT_VIEW) {
        if (strcmp(str_tolower(value),"Standard")==0)
          DEFAULT_VIEW=STANDARD;
        else if (strcmp(value,"Bestguess")==0)
          DEFAULT_VIEW=BEST_GUESS;
        else if (strcmp(value,"Userdefined")==0)
          DEFAULT_VIEW=USER_DEFINED;
        else
          ERROR;
      }
      
      else if (var==&CAMERA_X || var==&CAMERA_Y || var==&CAMERA_Z ||
               var==&LOOK_AT_X || var==&LOOK_AT_Y || var==&LOOK_AT_Z ||
               var==&SKY_X || var==&SKY_Y || var==&SKY_Z) {
        fVal=atof(value);
        *((F_TYPE *)var)=fVal;
      }

      else
        printf("Internal error at line %u! Contact the author.\n",lineNum);

    }

parsing_error:

    if (error) {
      printf(
        "Warning: parsing error in line %u of \"%s\". This line will be ignored.\n",
        lineNum,*fileName);
    }

  }
  fclose(config);
  printf("Read configuration file \"%s\".\n",*fileName);

/*
 * These are the text-based variables that didn't get initialized earlier
 * so that space wouldn't be wasted if they're specified in the config file.
 */

no_config:

  if (!PERIODIC_TABLE_FILE) New_String(&PERIODIC_TABLE_FILE,"periodic.tab");
  if (!FINISH) New_String(&FINISH,"Plastic");
  if (!FINISH_DECLARATION) New_String(&FINISH_DECLARATION,"\
  finish {\n\
    ambient 0.2\n\
    diffuse 0.7\n\
    reflection 0.0\n\
    brilliance 1.0\n\
    phong 0.3\n\
    phong_size 50\n\
    specular 0.0\n\
  }\n\n");

  if (POVRAY && (!POV_OPTIONS || !POV_WIDTH) ) {
    puts("Warning: to run POV-Ray from PovChem, you must define *all* the");
    puts("associated options in the configuration file!\n");
    POVRAY=NULL;
    RUN_POVRAY=FALSE;
  }

  if (!H_BOND_DONORS) New_String(&H_BOND_DONORS,"N,O");
  if (!H_BOND_ACCEPTORS) New_String(&H_BOND_ACCEPTORS,"N,O");
  
  return;
}


/*------------- set up Periodic Table --------------*/

void Parse_Periodic_Table(void)
{
  FILE *table;
  char line[256], tmpName[80], *tok;
  U_INT error=FALSE, lineNum=0;
  int Z;

  if ((table=fopen(PERIODIC_TABLE_FILE,"r"))==NULL) {
    printf("Can't read the periodic table from \"%s\"!\n",PERIODIC_TABLE_FILE);
    perror("Reason");
    exit(EXIT_FAILURE);
  }
  for (Z=0; Z<NUM_ELEMENTS; Z++) {
    PeriodicTable[Z].name=NULL;
    strcpy(PeriodicTable[Z].symbol,"");
    PeriodicTable[Z].isPresent=FALSE;
  }

  while (fgets(line,256,table)!=NULL) {
    lineNum++;
    error=FALSE;

    if (line[0]!='#' && !isBlank(line)) {

/*----- parse the line with really anal error checking ------*/

      Z=-1;
      tok=strtok(line," \t,");
      if (!tok) ERROR;
      Z=atoi(tok);
      tok=str_tolower(strtok(NULL," \t,"));
      if (!tok) ERROR;
      strcpy(tmpName,tok);
      if (Z==0 && strcmp(tmpName,"Unknown")!=0)
        ERROR;
      New_String(&(PeriodicTable[Z].name),tmpName);

      tok=str_tolower(strtok(NULL," \t,"));
      if (!tok) ERROR;
      strcpy(tmpName,tok);
      if (*tmpName=='\0') ERROR;
      strcpy(PeriodicTable[Z].symbol,tmpName);

      tok=strtok(NULL," \t,");
      if (!tok) ERROR;
      if ((PeriodicTable[Z].vdW=atof(tok))==0.0) ERROR;
        
      tok=strtok(NULL," \t,");
      if (!tok) ERROR;
      RED(Z)=atof(tok);
      if (RED(Z)>1.0 || RED(Z)<0.0) ERROR;

      tok=strtok(NULL," \t,");
      if (!tok) ERROR;
      GREEN(Z)=atof(tok);
      if (GREEN(Z)>1.0 || GREEN(Z)<0.0) ERROR;

      tok=strtok(NULL," \t,");
      if (!tok) ERROR;
      BLUE(Z)=atof(tok);
      if (BLUE(Z)>1.0 || BLUE(Z)<0.0) ERROR;
      
    }

parsing_error:
    if (error) {
      printf("Warning: error parsing line %u of \"%s\".\n",lineNum,
             PERIODIC_TABLE_FILE);
      if (Z!=-1) strcpy(PeriodicTable[Z].symbol,"?");
    }
  }

  fclose(table);
  printf("Read periodic table \"%s\".\n",PERIODIC_TABLE_FILE);

/*  for (Z=0;Z<NUM_ELEMENTS;Z++) {
    if (PeriodicTable[Z].name) {
      printf("%u %s %s %.2f %.2f %.2f %.2f\n",Z,PeriodicTable[Z].name,
             PeriodicTable[Z].symbol,PeriodicTable[Z].vdW,
             RED(Z),GREEN(Z),BLUE(Z));
    }
  }*/

  Make_H_Bonder_List(H_BOND_DONORS);
  Make_H_Bonder_List(H_BOND_ACCEPTORS);

  return;
}


/*-------- Search the PeriodicTable to determine Z from the element symbol --------*/

U_INT Determine_Z_From_PDB(char *line)
{
  char *pos=NULL, symbol[3];
  U_INT Z=0, i, try;

  for (try=0;try<2;try++) {
    if (try==0) {                 /* first try new PDB format field */
      if (strlen(line)>=78) {
        STRNCPY0(symbol,line+76,2);
        if (*symbol==' ')
          pos=symbol+1;
        else
          pos=symbol;
      } else
        try=1;
    }
    if (try==1) {                 /* if not, try to determine from atom name */
      pos=strpbrk(line+12,
                  "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz");
      if (pos<(line+16)) {
        STRNCPY0(symbol,pos,2);
        if (pos==line+15)
          symbol[1]='\0';
        pos=symbol;
      } else
        try=2;
    }
    if (try<2) {               /* try to determine Z from 2 or 1-char string */
      *pos=toupper(*pos);
      if (isalpha(*(pos+1))) {
        *(pos+1)=tolower(*(pos+1));
      } else 
        *(pos+1)='\0';
      for (i=1; i<=NUM_ELEMENTS; i++) {
        if (strncmp(pos,PeriodicTable[i].symbol,2)==0) {
          Z=i;
          try=2;
          break;
        }
      }
      if (Z==0) {
        for (i=1; i<=NUM_ELEMENTS; i++) {
          if (strncmp(pos,PeriodicTable[i].symbol,1)==0) {
            Z=i;
            try=2;
            break;
          }
        }
      }
    }
  }
  if (Z==0) {
    printf("WARNING: Can't determine element of atom serial #%5s\n",line+6);
    printf("Assigning vdW %.5g, color <%.5g, %.5g, %.5g>.\n",
           PeriodicTable[0].vdW,
           PeriodicTable[0].r,PeriodicTable[0].g,PeriodicTable[0].b);
  }
  return Z;
}

F_TYPE Determine_Radius(U_INT Z)
{
  switch (model) {
    case CPK:
      return PeriodicTable[Z].vdW;
      break;
    case BALLnSTICK:
      if (ballSizeModel==CONSTANT)
        return ballRadius;
      else
        return ballRadius*PeriodicTable[Z].vdW;
      break;
    case CYLINDER: default:
      return ballRadius;
  }
}

void Assign_Colors_To_Tags(U_INT *submodel)
{
  tagType *tag;
  U_INT i, nTags;
  F_TYPE deg;
  
  if (submodel==&sphereColorModel) {
    tag=firstSphereTag;
    nTags=nSphereTags;
  } else {
    tag=firstCylinderTag;
    nTags=nCylinderTags;
  }
  printf("Coloring %u different tags...\n\n",nTags);
  for (i=0;i<nTags;i++) {
    deg=(360.0/nTags)*i;
    if (deg>=0.0 && deg<120.0) {
      tag->r=1.0-deg/120.0;
      tag->g=1.0-tag->r;
    } else if (deg>=120.0 && deg<240.0) {
      tag->g=1.0-(deg-120.0)/120.0;
      tag->b=1.0-tag->g;
    } else {
      tag->b=1.0-(deg-240.0)/120.0;
      tag->r=1.0-tag->b;
    }
    tag=tag->next;
  }
  return;
}


/*---------- POV material (finish and pigment) statments and declarations ----------*/

void Write_Color_and_Atom_Declarations(void)
{
  U_INT Z, i, nTags;
  tagType *tag;

  if (nOrders[H_BOND]) {
    fputs("// Unfortunately, if the atoms involved in hydrogen-bonding are\n",pov);
    fputs("// resized, the hydrogen bonds will not be spaced properly. You\n",pov);
    fputs("// will have to re-run PovChem and change the sizes there.\n\n",pov);
  }

  if (model==CPK ||
      (model==BALLnSTICK && ballSizeModel==PROPORTIONAL) ) {
    if (model==BALLnSTICK)
      fprintf(pov,"#declare Proportionality_Constant = %.*g\n",SF,ballRadius);
    for (Z=0;Z<NUM_ELEMENTS;Z++) {
      if (PeriodicTable[Z].isPresent) {
        fprintf(pov,"#declare %s_Radius = %.*g",
                    PeriodicTable[Z].name,SF,PeriodicTable[Z].vdW);
        if (model==BALLnSTICK)
          fputs("*Proportionality_Constant\n",pov);
        else
          fputs("\n",pov);
      }
    }
  } else
    fprintf(pov,"#declare All_Atom_Radius = %.*g\n",SF,Determine_Radius(0));
  fputs("\n",pov);

  if (sphereColorModel==BY_ELEMENT || cylinderColorModel==BY_ELEMENT) {
    for (Z=0;Z<NUM_ELEMENTS;Z++) {
      if (PeriodicTable[Z].isPresent) {
        fprintf(pov,"#declare Atom_%s_Color = color rgb <%.5g, %.5g, %.5g>\n",
                    PeriodicTable[Z].name,RED(Z),GREEN(Z),BLUE(Z));
      }
    }
    fputs("\n",pov);
  }

  if (sphereColorModel==ALL_THE_SAME)
    fprintf(pov,"#declare All_Atom_Color = color rgb <%.5g, %.5g, %.5g>\n\n",
                allSphereColor.r,allSphereColor.g,allSphereColor.b);

  if (sphereColorModel==BY_PDB) {
    if (ballSizeModel==CONSTANT) {
      fputs("#declare Atom_All =\n  sphere {\n",pov);
      fputs("    <0,0,0>, All_Atom_Radius\n",pov);
      fprintf(pov,"    finish { %s }\n  }\n\n",FINISH);
    } else {
      for (Z=0;Z<NUM_ELEMENTS;Z++) {
        if (PeriodicTable[Z].isPresent) {
          fprintf(pov,"#declare Atom_%s =\n  sphere {\n",PeriodicTable[Z].name);
          fprintf(pov,"    <0,0,0>, %s_Radius\n",PeriodicTable[Z].name);
          fprintf(pov,"    finish { %s }\n  }\n",FINISH);
        }
      }
      fputs("\n",pov);
    }
    tag=firstSphereTag;
    nTags=nSphereTags;
    for (i=0;i<nTags;i++) {
      fprintf(pov,"#declare Atom_%s_Color = color rgb <%.5g, %.5g, %.5g>\n",
                  tag->name,tag->r,tag->g,tag->b);
      tag=tag->next;
    }

  } else if (model==CPK ||
             sphereColorModel==BY_ELEMENT || ballSizeModel==PROPORTIONAL) {
    for (Z=0;Z<NUM_ELEMENTS;Z++) {
      if (PeriodicTable[Z].isPresent) {
        fprintf(pov,"#declare Atom_%s =\n  sphere {\n",PeriodicTable[Z].name);

        if (model==CYLINDER ||
            (model==BALLnSTICK && ballSizeModel==CONSTANT) )
          fputs("    <0,0,0>, All_Atom_Radius\n",pov);
        else
          fprintf(pov,"   <0,0,0>, %s_Radius\n",PeriodicTable[Z].name);

        fputs(      "    texture {\n",pov);
        switch (sphereColorModel) {
          case BY_ELEMENT:
            fprintf(pov,"      pigment { color Atom_%s_Color }\n",
                        PeriodicTable[Z].name);
            break;
          case ALL_THE_SAME:
            fputs(      "      pigment { color All_Atom_Color }\n",pov);
            break;
        }
        fprintf(pov,"      finish { %s }\n    }\n  }\n",FINISH);
      }
    }

  } else if (sphereColorModel==ALL_THE_SAME) {
    fputs(      "#declare Atom_All =\n  sphere {\n",pov);
    fputs(      "    <0,0,0>, All_Atom_Radius\n",pov);
    fputs(      "    texture {\n",pov);
    fputs(      "      pigment { color All_Atom_Color }\n",pov);
    fprintf(pov,"      finish { %s }\n    }\n  }\n",FINISH);
  }

  fputs("\n",pov);
  return;
}

void Write_Bond_Declarations(void)
{
  U_INT i;
  tagType *tag;
  F_TYPE factors[5];

  factors[0]=0.0;
  factors[1]=1.0;
  factors[2]=DOUBLE_FACTOR;
  factors[3]=TRIPLE_FACTOR;
  factors[4]=HIGH_ORDER_FACTOR;

  if (nOrders[H_BOND]) {
    fprintf(pov,"#declare H_Bond_Radius = %.*g\n",SF,HBondRadius);
    fprintf(pov,"#declare H_Bond_Color = color rgb <%.5g, %.5g, %.5g>\n",
                HBondColor.r,HBondColor.g,HBondColor.b);
    fputs(      "#declare H_Bond =\n  cylinder {\n",pov);
    fputs(      "    <0,0,0>, <1,0,0>, H_Bond_Radius\n",pov);
    fputs(      "    texture {\n",pov);
    fputs(      "      pigment { H_Bond_Color }\n",pov);
    fprintf(pov,"      finish { %s }\n    }\n  }\n\n",FINISH);
  }

  if (model==CYLINDER) 
    fputs("#declare All_Bond_Radius = All_Atom_Radius\n\n",pov);
  else {
    if (cylinderSizeModel==CONSTANT)
      fprintf(pov,"#declare All_Bond_Radius = %.*g\n\n",SF,cylinderRadius);
    else {
      fprintf(pov,"#declare Single_Bond_Radius = %.*g\n",SF,cylinderRadius);
      for (i=DOUBLE;i<=HIGH_ORDER;i++) {
        if (nOrders[i])
          fprintf(pov,"#declare %sBond_Radius = Single_Bond_Radius*%.*g\n",
                     orderTypes[i],SF,factors[i]);
      }
      fputs("\n",pov);
    }
  }

  if (cylinderColorModel==ALL_THE_SAME)
    fprintf(pov,"#declare All_Bond_Color = color rgb <%.5g, %.5g, %.5g>\n\n",
                allCylinderColor.r,allCylinderColor.g,allCylinderColor.b);

  fputs(      "#declare Bond =\n  cylinder {\n",pov);
  
  if (cylinderSizeModel==CONSTANT)
    fputs(      "    <0,0,0>, <1,0,0>, All_Bond_Radius\n",pov);
  else
    fputs(      "    <0,0,0>, <1,0,0>, 1\t\t// leave as 1 - will be scaled later.\n",pov);
  
  if (cylinderColorModel==ALL_THE_SAME) {
    fputs(      "    texture {\n",pov);
    fprintf(pov,"      pigment { All_Bond_Color }\n");
    fprintf(pov,"      finish { %s }\n    }\n",FINISH);
  } else
    fprintf(pov,"    finish { %s }\n",FINISH);
    
  fputs(      "  }\n\n",pov);

  if (cylinderColorModel==BY_PDB && firstCylinderTag!=firstSphereTag) {
    tag=firstCylinderTag;
    for (i=0;i<nCylinderTags;i++) {
      fprintf(pov,"#declare Bond_%s_Color = color rgb <%.5g, %.5g, %.5g>\n",
                  tag->name,tag->r,tag->g,tag->b);
      tag=tag->next;
    }
  }

  return;
}


/*--------- math to determine distance between two points, and bond angle
            between three Atoms                                     -------*/

F_TYPE Distance(F_TYPE x1, F_TYPE y1, F_TYPE z1,
                F_TYPE x2, F_TYPE y2, F_TYPE z2)
{
  return SQRT((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
}

F_TYPE Angle(U_INT A, U_INT B, U_INT C)
{
  F_TYPE BAdist, BCdist, BAdotBC;     /* ... determine A-B-C angle        */
                                      /* from cos(angle)=(BA.BC)/|BA||BC| */

  BAdist=Distance(X(A),Y(A),Z(A),X(B),Y(B),Z(B));
  BCdist=Distance(X(C),Y(C),Z(C),X(B),Y(B),Z(B));
  if (BAdist==0.0 || BCdist==0.0) return 0.0;
  BAdotBC=(X(A)-X(B))*(X(C)-X(B))+
          (Y(A)-Y(B))*(Y(C)-Y(B))+
          (Z(A)-Z(B))*(Z(C)-Z(B));
  return DEGREES(ACOS(BAdotBC/(BAdist*BCdist)));
}


/*------------- reads a (possibly null) string from stdin ------------*/

char *Get_User_Input(char *string)
{
  int i=0;
  
  do {
    scanf("%c",&string[i]);
  } while (string[i++]!='\n');
  string[i-1]='\0';
  return string;
}


/*------------ write the POV sphere with appropriate element & radius ----------*/

#define RAND_360 ((360.0/MAXINT)*random())

void Write_Atom(F_TYPE x, F_TYPE y, F_TYPE z,
                U_INT Z, char *tag)
{

  aveX+=x;            /* need these to set up camera location later... */
  aveY+=y;
  aveZ+=z;
  if (x<minX) minX=x;
  if (x>maxX) maxX=x;
  if (y<minY) minY=y;
  if (y>maxY) maxY=y;
  if (z<minZ) minZ=z;
  if (z>maxZ) maxZ=z;
  
  if (sphereColorModel==BY_PDB) {

    if (model==CPK ||
        (model==BALLnSTICK && ballSizeModel==PROPORTIONAL) )
      fprintf(inc,"object { Atom_%s\n",PeriodicTable[Z].name);
    else
      fputs(      "object { Atom_All\n",inc);
    fprintf(inc,"  pigment { color Atom_%s_Color }\n  ",tag);

  } else
    fprintf(inc,"object { Atom_%s ",tag);

/*  fprintf(inc,"\n  rotate <%.1f, %.1f, %.1f>\n  ",
               RAND_360,RAND_360,RAND_360); */
  fprintf(inc,"translate <%.*g, %.*g, %.*g> }\n",SF,x,SF,y,SF,z);
  if (model==CPK || sphereColorModel==BY_ELEMENT ||
      cylinderColorModel==BY_ELEMENT || ballSizeModel==PROPORTIONAL)
    PeriodicTable[Z].isPresent=TRUE;
  return;
}

void Parse_PDB_Line(char *line, PDBatomInfoType *atom, U_INT searchFlags)
{
   char temp[9];
   int i;

   if (searchFlags & FIND_SERIAL) {
     STRNCPY0(temp,line+6,5);
     atom->serial=atoi(temp);
   }

   if (searchFlags & FIND_Z)
     atom->Z=Determine_Z_From_PDB(line);

   if (searchFlags & FIND_RESNAME) {
     STRNCPY0(temp,line+17,3);
     for (i=0;temp[i]==' ';i++) {}
     strcpy(atom->resName,&temp[i]);
     for (i=2;(atom->resName[i]==' ' && i>=0);i--)
       atom->resName[i]='\0';
   }

   if (searchFlags & FIND_CHAINID) {
     atom->chainID=line[21];
     if (atom->chainID==' ') atom->chainID='_';
   }

   if (searchFlags & FIND_RESSEQ) {
     STRNCPY0(temp,line+22,4);
     atom->resSeq=atoi(temp);
   }

   if (searchFlags & FIND_COORDS) {
     STRNCPY0(temp,line+30,8);
     atom->x=(F_TYPE)atof(temp);
     STRNCPY0(temp,line+38,8);
     atom->y=(F_TYPE)atof(temp);
     STRNCPY0(temp,line+46,8);
     atom->z=(F_TYPE)(-atof(temp));  /* NOTE: POV Z axis is mirror of PDB's */
   }

   return;
}

char *Make_Tag_From_PDB_Info(PDBatomInfoType info, U_INT *submodel)
{
  static char tag[80], tmp[80];
  char defaultFields;
  U_INT length, *choice;

  if (submodel==&sphereColorModel) {
    choice=&sphereTagType;
    defaultFields=DEFAULT_PDB_FIELDS;
  } else {
    choice=&stickTagType;
    defaultFields=DEFAULT_STICK_PDB_FIELDS;
  }
  if (*choice==0) {
    if (isInteractive) {
      switch (model) {
        case CPK:
          strcpy(tmp,   "Atoms");
          strcpy(tmp+40,"     ");
          break;
        case CYLINDER:
          strcpy(tmp,   "Cylinders");
          strcpy(tmp+40,"         ");
          break;
        case BALLnSTICK: default:
          if (submodel==&sphereColorModel) {
            strcpy(tmp,   "Balls");
            strcpy(tmp+40,"     ");
          } else {
            strcpy(tmp,   "Sticks");
            strcpy(tmp+40,"      ");
          }
          break;
      }
      puts("");
      printf("Tag %s by: 1) Residue\n",tmp);
      printf("    %s     2) Chain\n",tmp+40);
      printf("    %s     3) Residue + Chain + Sequence number\n",tmp+40);
      printf("    %s     4) Residue + Element\n",tmp+40);
      printf("Choice [%u] ? ",defaultFields);
      *choice=(Get_User_Input(tmp)[0]=='\0' ? defaultFields : (char)atoi(tmp));
      if (*choice<1 || *choice>4) *choice=defaultFields;
      puts("");
    } else
      *choice=defaultFields;
  }
  switch (*choice) {
    case 1:
      strcpy(tag,info.resName);
      break;
    case 2:
      tag[0]=info.chainID;
      tag[1]='\0';
      break;
    case 3:
      strcpy(tag,info.resName);
      strcat(tag,"_");
      length=strlen(tag);
      tag[length]=info.chainID;
      tag[length+1]='_';
      sprintf(&(tag[length+2]),"%u",info.resSeq);
      break;
    case 4:
      strcpy(tag,PeriodicTable[info.Z].name);
      strcat(tag,"_");
      strcat(tag,info.resName);
      break;
  }
  return tag;
}

void Make_Tag_List(tagType **firstTag, U_INT *nTags)
{
  tagType *currentTag, *prevTag;
  char tmp[81], currentName[80];
  PDBatomInfoType atom;
  char isInList;

  for ( ; (fgets(tmp,81,pdb))!=NULL && strncmp(tmp,"END   ",6)!=0 ; ) {
    if (strncmp(tmp,"ATOM  ",6)==0 || strncmp(tmp,"HETATM",6)==0) {
      
      Parse_PDB_Line(tmp,&atom,
                     (FIND_Z | FIND_RESNAME | FIND_CHAINID | FIND_RESSEQ));
      if (nTags==&nSphereTags)
        strcpy(currentName,Make_Tag_From_PDB_Info(atom,&sphereColorModel));
      else
        strcpy(currentName,Make_Tag_From_PDB_Info(atom,&cylinderColorModel));
      isInList=FALSE;
      currentTag=*firstTag;
      prevTag=NULL;
      do {
        if (currentTag!=NULL) {
          if (strcmp(currentName,currentTag->name)==0) {
            isInList=TRUE;
          } else {
            prevTag=currentTag;
            currentTag=currentTag->next;
          }
        }
        if (currentTag==NULL) {
          if ((currentTag=(tagType *)malloc(sizeof(tagType)))==NULL) {
            puts("FATAL: Not enough memory for even one more tag!");
            exit(EXIT_FAILURE);
          } else {
            memAlloced+=sizeof(tagType);
            (*nTags)++;
            New_String(&(currentTag->name),currentName);
            currentTag->next=NULL;
            currentTag->r=currentTag->g=currentTag->b=0.0;
            isInList=TRUE;
            if (*firstTag==NULL) {
              printf("\nFound tags: %s",currentTag->name);
              *firstTag=currentTag;
            } else {
              prevTag->next=currentTag;
              printf(", %s",currentTag->name);
            }
          }
        }
      } while (isInList==FALSE);
    }
  }
  puts("\n");
  rewind(pdb);
  return;
}


/*-------------  to write out a POV cylinder for the bond -----------------*/

/*
Okay, I know this is a roundabout way of doing it, but it has to be this way
in order to be able to put any arbitrary object in place of a bond, just in
case I wanted to do that someday... :) Also, this allows split-color bonds
to be made out of a single cylinder with an appropriate color map.

Anyway, this assumes the bond is originally a cylinder from <0,0,0> to <1,0,0>
of proper radius. Then 1) color and finish it if necessary, 2) scale it up to
the right length, 3) rotate around two angles about the Y and Z axes to align
the bond along the proper vector, and 4) translate the thing into place.
*/

void Write_Cylinder(F_TYPE x1, F_TYPE y1, F_TYPE z1,
                    F_TYPE x2, F_TYPE y2, F_TYPE z2,
                    char *radiusType,
                    char *tag1, char *tag2,
                    F_TYPE vdW1, F_TYPE vdW2)
{
  F_TYPE length, dy, theta, phi, proportion;
  char type[10], tmp[80];
  
  length=Distance(x1,y1,z1,x2,y2,z2);
  if (length==0.0) return;

  if (radiusType==orderTypes[H_BOND])
    fputs(      "object { H_Bond\n",inc);
  else {
    fputs(      "object { Bond\n",inc);
    if (model==BALLnSTICK && cylinderSizeModel==PROPORTIONAL)
      fprintf(inc,"  scale <1, %sBond_Radius, %sBond_Radius>\n",
                  radiusType,radiusType);
  }

  if (radiusType==orderTypes[H_BOND])
    type[0]='\0';
  else if (cylinderColorModel==BY_PDB && firstCylinderTag!=firstSphereTag)
    strcpy(type,"Bond_");
  else
    strcpy(type,"Atom_");

  if (radiusType!=orderTypes[H_BOND] && cylinderColorModel!=ALL_THE_SAME) {
    if (strcmp(tag1,tag2)==0) {
      fprintf(inc,"  pigment { color %s%s_Color }\n",type,tag1);
    } else {
      fputs(      "  pigment {\n",inc);
      fputs(      "    gradient x\n",inc);
      fputs(      "    color_map {\n",inc);
      
      switch (stickSplit) {
      case SPLIT_EQUAL:
        fprintf(inc,"      [0.5 color %s%s_Color ]\n",type,tag1);
        fprintf(inc,"      [0.5 color %s%s_Color ]\n    }\n  }\n",type,tag2);
        break;
      case SPLIT_PROPORTIONAL:
        proportion=POW(vdW1,POWER)/
                (POW(vdW1,POWER)+POW(vdW2,POWER));
        fprintf(inc,"      [%.5g color %s%s_Color ]\n",
                    proportion,type,tag1);
        fprintf(inc,"      [%.5g color %s%s_Color ]\n    }\n  }\n",
                    proportion,type,tag2);
        break;
      case SPLIT_SHADED:
        fprintf(inc,"      [0.15 color %s%s_Color ]\n",type,tag1);
        fprintf(inc,"      [0.85 color %s%s_Color ]\n    }\n  }\n",type,tag2);
        break;
      }
    }
  }
  
  fprintf(inc,"  scale <%.*g, 1, 1>\n",SF,length);
  
  phi=ACOS((y2-y1)/length);
  sprintf(tmp,"  rotate <0, 0, %.*g>\n",SF,DEGREES(-phi)+90.0);
  if (strcmp(tmp,"  rotate <0, 0, 0>\n")!=0)
    fprintf(inc,"%s",tmp);

  dy=SQRT((x2-x1)*(x2-x1)+(z2-z1)*(z2-z1));
  if (dy!=0.0) {
    theta=ACOS((x2-x1)/dy);
    sprintf(tmp,"  rotate <0, %.*g, 0>\n",SF,
                (z2-z1)>0.0 ? DEGREES(-theta) : DEGREES(theta));
    if (strcmp(tmp,"  rotate <0, 0, 0>\n")!=0)
      fprintf(inc,"%s",tmp);
  }

  fprintf(inc,"  translate <%.*g, %.*g, %.*g>  }\n",SF,x1,SF,y1,SF,z1);
  return;
}


/*---------- find indexes to the "Atoms" array of the atoms given
             by their serial number ----------------------------*/

#define FIRST 1
#define SECOND 2

char Find_Indexes(U_INT serNumA, U_INT serNumB,
                          U_INT *indexA, U_INT *indexB)
{
  U_INT i, found=0;
  
  if (serNumA>0 && Atoms[serNumA-1].serNum==serNumA) {
    *indexA=serNumA-1;
    found|=FIRST;                      /* serial #'s usually in order. */
  }
  if (serNumB>0 && Atoms[serNumB-1].serNum==serNumB) {
    *indexB=serNumB-1;
    found|=SECOND;
  }
  for ( i=0; (i<nAtoms && found<(FIRST|SECOND)); i++ ) {
    if (!(found&FIRST) && Atoms[i].serNum==serNumA) {
      *indexA=i;
      found|=FIRST;                   /* but if not, search the long way... */
    }
    else if (!(found&SECOND) && Atoms[i].serNum==serNumB) {
      *indexB=i;
      found|=SECOND;
    }
  }
  if (found!=(FIRST|SECOND)) {
    printf("Unable to find both atoms %u and %u!\n",serNumA,serNumB);
    return 1;
  }
  return 0;
}


/* -------------- update the bond order in memory ----------------------*/

void Add_Bond(U_INT serNumA, U_INT serNumB, unsigned char order)
{
  U_INT indexA, indexB, a, b, *tmpA=0, *tmpB=0;
  unsigned char *tmpO;
  
  if (Find_Indexes(serNumA,serNumB,&indexA,&indexB)) {
    printf("Bond between serial #'s %u and %u ignored.\n",serNumA,serNumB);
  } else {
    for (a=0; a<MAX_BONDS_PER_ATOM &&
              (*(tmpA=&BoundTo(indexA,a)))!=indexB &&
              *tmpA!=nAtoms    ; a++) {}
    if (a==MAX_BONDS_PER_ATOM) {
      printf("\nFatal! MAX_BONDS_PER_ATOM exceeded for atom #%u.\n",indexA);
      BOMB;
    }
    for (b=0; b<MAX_BONDS_PER_ATOM &&
              (*(tmpB=&BoundTo(indexB,b)))!=indexA &&
              *tmpB!=nAtoms    ; b++) {}
    if (b==MAX_BONDS_PER_ATOM) {
      printf("\nFatal! MAX_BONDS_PER_ATOM exceeded for atom #%u.\n",indexB);
      BOMB;
    }
    *tmpA=indexB;
    *tmpB=indexA;
    if (order==H_BOND) {
      BondOrder(indexA,a)=BondOrder(indexB,b)=H_BOND;
    } else {
      tmpO=&(BondOrder(indexA,a));
      if ( ++(*tmpO) > HIGH_ORDER)
        *tmpO=HIGH_ORDER;
      tmpO=&(BondOrder(indexB,b));
      if ( ++(*tmpO) > HIGH_ORDER)
        *tmpO=HIGH_ORDER;
    }
  }
  return;
}


/*-- function to find the positions of ALL bond ends and write cylinders ---*/

void Write_Bonds(void)
{
  char *name1, *name2, *radiusType;
  unsigned char order;
  U_INT index1, index2, i, j, nDashes;
  F_TYPE proportion, proportion2, atomRadius1, atomRadius2;
  F_TYPE a1, b1, c1, a2, b2, c2, length, dashSep;

  nOrders[SINGLE]=nOrders[DOUBLE]=nOrders[TRIPLE]=
   nOrders[HIGH_ORDER]=nOrders[H_BOND]=0;

  for (index1=0; index1<nAtoms; index1++) {
  for (j=0; j<MAX_BONDS_PER_ATOM &&
            (index2=BoundTo(index1,j))!=nAtoms; j++) {
  if (index2>index1) {

  order=BondOrder(index1,j);
  (nOrders[order])++;
  radiusType=orderTypes[order];

  if (order==H_BOND) {

    atomRadius1=Determine_Radius(Atoms[index1].Z);
    atomRadius2=Determine_Radius(Atoms[index2].Z);
    length=Distance(X(index1),Y(index1),Z(index1),
                    X(index2),Y(index2),Z(index2));
    proportion=(atomRadius1)/length;
    proportion2=(atomRadius2)/length;
    a1=X(index1)+proportion*(X(index2)-X(index1));
    b1=Y(index1)+proportion*(Y(index2)-Y(index1));
    c1=Z(index1)+proportion*(Z(index2)-Z(index1));
    a2=X(index2)+proportion2*(X(index1)-X(index2));
    b2=Y(index2)+proportion2*(Y(index1)-Y(index2));
    c2=Z(index2)+proportion2*(Z(index1)-Z(index2));
    dashSep=1.9/(2.0*HBondDashes+1.0);
    nDashes=Round((length/dashSep-1.0)/2.0);
    for (i=0;i<nDashes;i++) {
      Write_Cylinder(a1+(a2-a1)*(((F_TYPE)(2*i+1))/(2*nDashes+1)),
                     b1+(b2-b1)*(((F_TYPE)(2*i+1))/(2*nDashes+1)),
                     c1+(c2-c1)*(((F_TYPE)(2*i+1))/(2*nDashes+1)),
                     a1+(a2-a1)*(((F_TYPE)(2*i+2))/(2*nDashes+1)),
                     b1+(b2-b1)*(((F_TYPE)(2*i+2))/(2*nDashes+1)),
                     c1+(c2-c1)*(((F_TYPE)(2*i+2))/(2*nDashes+1)),
                     radiusType,"H_Bond","H_Bond",0,0);
    }
    
  } else {
  
    switch (cylinderColorModel) {
      case ALL_THE_SAME:
        name1=name2="All";
        break;
      case BY_ELEMENT: case BY_PDB: default:
        if (cylinderColorModel==BY_PDB) {
          if (model==BALLnSTICK) {
            name1=cylinderEndTags[index1];
            name2=cylinderEndTags[index2];
          } else {
            name1=TAG(index1);
            name2=TAG(index2);
          }
        } else {
          name1=NAME(index1);
          name2=NAME(index2);
        }
        break;
    }
    Write_Cylinder(X(index1),Y(index1),Z(index1),
                   X(index2),Y(index2),Z(index2),
                   radiusType,name1,name2,VDW(index1),VDW(index2));
  }
  
  }}}
  
  printf(
    "Wrote %u single, %u double, %u triple, %u higher order, and %u hydrogen bonds.\n",
    nOrders[SINGLE],nOrders[DOUBLE],nOrders[TRIPLE],
    nOrders[HIGH_ORDER],nOrders[H_BOND]);
  return;
}


/*------------------- Hydrogen bond functions ----------------------*/

#define HYDROGEN 1
#define NITROGEN 7
#define OXYGEN 8

void Find_Hydrogen_Bonds(void)       /* H-Bond: D-H...A-R */
{
  U_INT D, H, A, R, DZ, AZ, indexH;
  F_TYPE HAdist, DHAangle;
  unsigned char tmp;
  
  for (D=0;D<nAtoms;D++) {         /* find all H-bond donors... */
   for (DZ=0; DZ<nDonors && HBondDonors[DZ]!=ZNUM(D) ; DZ++) {}
   if (DZ<nDonors) {              /* if D is in donor list */
    
    for (H=0;H<MAX_BONDS_PER_ATOM;H++) {  /* find any H's attached to them... */
     if ( (tmp=BondOrder(D,H))>ZERO &&
          tmp<H_BOND &&         /* (by covalent bond) */
          ZNUM(indexH=BoundTo(D,H))==HYDROGEN ) {
        
      for (A=0;A<nAtoms;A++) {  /* find all relevant H-bond acceptors */
       for (AZ=0; AZ<nAcceptors && HBondAcceptors[AZ]!=ZNUM(A) ; AZ++) {}
       if ((A!=D) && AZ<nAcceptors) {       /* if A is in acceptor list */
            
        HAdist=Distance(X(A),Y(A),Z(A),X(indexH),Y(indexH),Z(indexH));
        if (HAdist>=HBondMinDist && HAdist<=HBondMaxDist) {  /* check distance */
        
         DHAangle=Angle(D,indexH,A);
         if (DHAangle>=HBondMinDHA) {           /* check DHA angle */
                
          for (R=0;R<MAX_BONDS_PER_ATOM;R++) { /* check angle to any
                                                  covalent bonds to Acceptor  */
                                                            
           if ( (tmp=BondOrder(A,R))>ZERO &&     /* if too small, no go */
                tmp<H_BOND &&
                Angle(indexH,A,BoundTo(A,R))<HBondMinHAR)
             R=MAX_BONDS_PER_ATOM+1;
          }
          if (R==MAX_BONDS_PER_ATOM) {
           Add_Bond(SERNUM(indexH),SERNUM(A),H_BOND);
          }
         }
        }
       }
      }
     }
    }
   }
  }
  return;
}


/*------------- deals with user-input color values -----------*/

void Parse_Colors(char *input, Color *color)
{
  char okay, *entry, *endptr;
  Color defCol;
  
  defCol.r=color->r;
  defCol.g=color->g;
  defCol.b=color->b; 
  do {
    if (*input=='\0') {
      color->r=defCol.r;
      color->g=defCol.g;
      color->b=defCol.b;
      return;
    }
    okay=TRUE;

    entry=strtok(input," ,");
    if (entry) {
      color->r=(F_TYPE)strtod(entry,&endptr);
      if (endptr==(char *)entry || color->r<0.0 || color->r>1.0) {
        okay=FALSE;
        goto problem;
      }
    } else {
      okay=FALSE;
      goto problem;
    }

    entry=strtok(NULL," ,");
    if (entry) {
      color->g=(F_TYPE)strtod(entry,&endptr);
      if (endptr==(char *)entry || color->g<0.0 || color->g>1.0) {
        okay=FALSE;
        goto problem;
      }
    } else {
      okay=FALSE;
      goto problem;
    }

    entry=strtok(NULL," ,");
    if (entry) {
      color->b=(F_TYPE)strtod(entry,&endptr);
      if (endptr==(char *)entry || color->b<0.0 || color->b>1.0) {
        okay=FALSE;
      }
    } else {
      okay=FALSE;
    }

problem:
    if (!okay) {
      puts("You must give me a comma or space-separated list of three");
      puts("values between 0.0 and 1.0.\n");
      printf("What color [%.5g, %.5g, %.5g] ? ",defCol.r,defCol.g,
             defCol.b);
      Get_User_Input(input);    
    }
  } while (okay==FALSE);
  return;
}


/*
 * This function rotates a point around an arbitrary axis through the origin.
 */

void Apply_Rotation(F_TYPE axisX, F_TYPE axisY, F_TYPE axisZ,
                    F_TYPE radians,
                    F_TYPE *pointX, F_TYPE *pointY, F_TYPE *pointZ)
{
  F_TYPE lengthP, lengthA, scale, lengthPC;
  Vector C, H, V;
  
/*
 * C is a point along A (the axis). The plane perpendicular to A that
 * intersects axis A at point C also contains P, the point to be rotated.
 * Unit vector H points from C to P, and defines the "horizontal" axis in the
 * plane. V, the vertical unit length in the plane, is the cross-product of H
 * and A. R, the point resulting from rotating P theta radians around A, is
 * then C + |PC|*V*sin(theta) + |PC|*H*cos(theta).
 */
 
  lengthA=Distance(0.0,0.0,0.0,axisX,axisY,axisZ);
  lengthP=Distance(0.0,0.0,0.0,*pointX,*pointY,*pointZ);
  if (lengthP==0.0) return;
  scale=((*pointX)*axisX+(*pointY)*axisY+(*pointZ)*axisZ)/(lengthA*lengthA);
  C.x=axisX*scale;
  C.y=axisY*scale;
  C.z=axisZ*scale;
  H.x=(*pointX)-C.x;
  H.y=(*pointY)-C.y;
  H.z=(*pointZ)-C.z;
  lengthPC=Distance(C.x,C.y,C.z,*pointX,*pointY,*pointZ);
  if (lengthPC==0.0) { return; }
  V.x=axisY*H.z-axisZ*H.y;
  V.y=axisZ*H.x-axisX*H.z;
  V.z=axisX*H.y-axisY*H.x;
  H.x/=lengthPC;
  H.y/=lengthPC; /* make H a unit vector */
  H.z/=lengthPC;
  scale=Distance(0.0,0.0,0.0,V.x,V.y,V.z);
  V.x/=scale;
  V.y/=scale;    /* make V a unit vector */
  V.z/=scale;
  *pointX=C.x+V.x*lengthPC*SIN(radians)+H.x*lengthPC*COS(radians);
  *pointY=C.y+V.y*lengthPC*SIN(radians)+H.y*lengthPC*COS(radians);
  *pointZ=C.z+V.z*lengthPC*SIN(radians)+H.z*lengthPC*COS(radians);
  return;
}


/*---------------- camera set up routines --------------------*/

#define MULT (2.0)   /* to guess camera distance, defaults to  */
                     /* MULT * maximum spread in X, Y, or Z    */

void Set_Up_View(void)
{
  F_TYPE Xspread, Yspread, Zspread, length, theta, lengthS;
  char user[80], defaultBG, choice;
  Vector lineOfSight;
  
  Xspread=maxX-minX;
  Yspread=maxY-minY;
  Zspread=maxZ-minZ;

  if (isInteractive) {
    puts("\nDo you want: 1) Best-guess View");
    puts(  "             2) Standard View along -Z, +Y up, +X right");
    puts(  "             3) User-defined View");
    printf("Choice [%u] ? ",DEFAULT_VIEW);
    viewType=( (Get_User_Input(user)[0]=='\0') ? DEFAULT_VIEW : atoi(user) );
    if (viewType<1 || viewType>3) viewType=DEFAULT_VIEW;
    puts("");
  } else
    viewType=DEFAULT_VIEW;

  if (viewType!=USER_DEFINED) {
    if (isInteractive) {
      puts(  "Make the image orientation: 1) Landscape");
      puts(  "                            2) Portrait");
      printf("Choice [%u] ? ",DEFAULT_ORIENTATION);
      imageOrientation=
        ((Get_User_Input(user)[0]=='\0') ? DEFAULT_ORIENTATION : atoi(user));
      if (imageOrientation<1 || imageOrientation>2)
        imageOrientation=DEFAULT_ORIENTATION;
      puts("");
    } else
      imageOrientation=DEFAULT_ORIENTATION;
  }
  
  if (viewType==USER_DEFINED) {
    if (WINDOW_UP==0.0 || WINDOW_RIGHT==0.0) {
      puts("Fatal: Window has zero width or height! Aborting.\n");
      exit(EXIT_FAILURE);
    }
  } else {
    CAMERA_X=LOOK_AT_X=aveX;       /* camera will be moved later */
    CAMERA_Y=LOOK_AT_Y=aveY;
    CAMERA_Z=LOOK_AT_Z=aveZ;
    DIRECTION=2.0;
    SKY_X=SKY_Y=SKY_Z=0.0;
    unitRight.x=unitRight.y=unitRight.z=0.0;
    switch (imageOrientation) {
      case LANDSCAPE:
        WINDOW_UP=1.0;
        WINDOW_RIGHT=1.333;
        break;
      case PORTRAIT:
        WINDOW_UP=1.333;
        WINDOW_RIGHT=1.0;
        break;
    }
  }

  switch (viewType) {
    case BEST_GUESS:
      if (Yspread>Xspread) {
        if (Zspread>Yspread) {                 /* Z>Y>X */
          if (imageOrientation==LANDSCAPE) {
            CAMERA_X+=Zspread*MULT;
            SKY_Y=1.0;
            unitRight.z=1.0;
          } else {
            CAMERA_X+=Zspread*MULT;
            SKY_Z=1.0;
            unitRight.y=-1.0;
          }
        } else {
          if (Zspread>Xspread) {               /* Y>Z>X */
            if (imageOrientation==LANDSCAPE) {
              CAMERA_X+=Yspread*MULT;
              SKY_Z=1.0;
              unitRight.y=-1.0;
            } else {
              CAMERA_X+=Yspread*MULT;
              SKY_Y=1.0;
              unitRight.z=1.0;
            }
          } else {                             /* Y>X>Z */
            if (imageOrientation==LANDSCAPE) {
              CAMERA_Z-=Yspread*MULT;
              SKY_X=1.0;
              unitRight.y=-1.0;
            } else {
              CAMERA_Z-=Yspread*MULT;
              SKY_Y=1.0;
              unitRight.x=1.0;
            }
          }
        }
      } else {
        if (Zspread>Xspread) {                 /* Z>X>Y */
          if (imageOrientation==LANDSCAPE) {
            CAMERA_Y+=Zspread*MULT;
            SKY_X=1.0;
            unitRight.z=-1.0;
          } else {
            CAMERA_Y+=Zspread*MULT;
            SKY_Z=1.0;
            unitRight.x=1.0;
          }
        } else {
          if (Zspread>Yspread) {               /* X>Z>Y */
            if (imageOrientation==LANDSCAPE) {
              CAMERA_Y+=Xspread*MULT;
              SKY_Z=1.0;
              unitRight.x=1.0;
            } else {
              CAMERA_Y+=Xspread*MULT;
              SKY_X=1.0;
              unitRight.z=-1.0;
            }
          } else {                             /* X>Y>Z */
            if (imageOrientation==LANDSCAPE) {
              CAMERA_Z-=Xspread*MULT;
              SKY_Y=1.0;
              unitRight.x=1.0;
            } else {
              CAMERA_Z-=Xspread*MULT;
              SKY_X=1.0;
              unitRight.y=-1.0;
            }
          }
        }
      }
      break;
    
    case STANDARD:
      if (imageOrientation==LANDSCAPE)
        CAMERA_Z-=Xspread*MULT;
      else
        CAMERA_Z-=Yspread*MULT;
      SKY_Y=1.0;                          /* look along Z, sky along y */
      unitRight.x=1.0;
      break;
    
    case USER_DEFINED:
      lineOfSight.x=LOOK_AT_X-CAMERA_X;
      lineOfSight.y=LOOK_AT_Y-CAMERA_Y;
      lineOfSight.z=LOOK_AT_Z-CAMERA_Z;
      unitRight.x=SKY_Y*lineOfSight.z-SKY_Z*lineOfSight.y;
      unitRight.y=SKY_Z*lineOfSight.x-SKY_X*lineOfSight.z;
      unitRight.z=SKY_X*lineOfSight.y-SKY_Y*lineOfSight.x;
      length=Distance(0.0,0.0,0.0,unitRight.x,unitRight.y,unitRight.z);
      if (length==0.0) {
        puts("Fatal: error due to contradictory or meaningless camera,");
        puts("look at, or sky vectors. Aborting.\n");
        exit(EXIT_FAILURE);
      }
      unitRight.x/=length;
      unitRight.y/=length;
      unitRight.z/=length;

      length=Distance(0.0,0.0,0.0,lineOfSight.x,lineOfSight.y,lineOfSight.z);
      lengthS=Distance(0.0,0.0,0.0,SKY_X,SKY_Y,SKY_Z);
      SKY_X/=lengthS;
      SKY_Y/=lengthS;  /* make Sky unit length */
      SKY_Z/=lengthS;

      theta=ACOS(
        (SKY_X*lineOfSight.x+SKY_Y*lineOfSight.y+SKY_Z*lineOfSight.z)/
        (length) );
      if ( DEGREES(theta)>90.0001 || DEGREES(theta)<89.9999 ) {
/*
        puts("\nWarning: sky vector not perpendicular to line of sight!");
        puts("Applying correction.\n");
*/
        Apply_Rotation(unitRight.x,unitRight.y,unitRight.z,
                       theta-M_PI_2,
                       &SKY_X,&SKY_Y,&SKY_Z);
      }

      break;
  }
  
  if (isInteractive) {
    puts(  "Do you want the background: 1) Black");
    puts(  "                            2) White");
    puts(  "                            3) Dark blue");
    puts(  "                            4) Other");
    if (DEFAULT_BACKGROUND_RED+
        DEFAULT_BACKGROUND_GREEN+
        DEFAULT_BACKGROUND_BLUE==0.0)
      defaultBG=1;
    else
      defaultBG=4;
    printf("Choice [%u] ? ",defaultBG);
    choice=
      ( (Get_User_Input(user)[0]=='\0') ? defaultBG : atoi(user) );
    if (choice<1 || choice>4) choice=defaultBG;
    puts("");
  } else
    choice=4;
  switch (choice) {
    case 1:
      Background.r=Background.g=Background.b=0.0;
      break;
    case 2:
      Background.r=Background.g=Background.b=1.0;
      break;
    case 3:
      Background.r=0.0;
      Background.g=0.0;
      Background.b=0.19608;
      break;
    case 4: default:
      Background.r=DEFAULT_BACKGROUND_RED;
      Background.g=DEFAULT_BACKGROUND_GREEN;
      Background.b=DEFAULT_BACKGROUND_BLUE;
      if (isInteractive) {
        printf("What color [%.5g, %.5g, %.5g] ? ",DEFAULT_BACKGROUND_RED,
               DEFAULT_BACKGROUND_GREEN,DEFAULT_BACKGROUND_BLUE);
          Parse_Colors(Get_User_Input(user),&Background);
          puts("");
        }
      break;
  }
  
  return;
}


/*------------- writes the main header .pov file ------------------*/

void Write_POV_Header(U_INT whichEye, char *fileName)
{
  F_TYPE cameraDist, eyeDist, XZproj, rot1_Z, rot2_X, rot3_Y;
  Vector lineOfSight, eye, eyeOffset, invSky;
  char tmp[256];

  if ( (pov=fopen(fileName,"w"))==0 ) {
    printf("Can't write pov header file \"%s\"\n",fileName);
    perror("Reason");
    return;
  }
  fprintf(pov,"// This file was created by PovChem version %s\n\n",VERSION);

/*
 * If viewType is USER_DEFINED, then the camera will be initially at the
 * origin, then rotated, then translated into place. Otherwise the camera is
 * put directly into place.
 *
 * This function calculates all appropriate camera positions and angles,
 * and light positions, including slight camera offset for stereo images.
 */
 
  cameraDist=Distance(CAMERA_X,CAMERA_Y,CAMERA_Z,LOOK_AT_X,LOOK_AT_Y,LOOK_AT_Z);
  eye.x=CAMERA_X;
  eye.y=CAMERA_Y;
  eye.z=CAMERA_Z;

  if (whichEye!=MONO) {
    eyeOffset.x=unitRight.x*cameraDist*(STEREO_SEPARATION_FACTOR/2.0);
    eyeOffset.y=unitRight.y*cameraDist*(STEREO_SEPARATION_FACTOR/2.0);
    eyeOffset.z=unitRight.z*cameraDist*(STEREO_SEPARATION_FACTOR/2.0);
    if (whichEye==LEFT_EYE) {
      eyeOffset.x=-eyeOffset.x;
      eyeOffset.y=-eyeOffset.y;
      eyeOffset.z=-eyeOffset.z;
    }
    eye.x+=eyeOffset.x;
    eye.y+=eyeOffset.y;       /* change view position if for stereo image */
    eye.z+=eyeOffset.z;
  }
  eyeDist=Distance(eye.x,eye.y,eye.z, LOOK_AT_X,LOOK_AT_Y,LOOK_AT_Z);


/*
 * Always set up the basic camera view window first, with camera placed at
 * the origin, looking along +Z, sky along +Y.
 */
 
  fprintf(pov,"#declare Camera_Position = < %.6g, %.6g, %.6g >\n\n",
              eye.x,eye.y,eye.z);

  fputs(      "camera{\n",pov);
  fputs(      "  location < 0, 0, 0 >\n",pov);
  fprintf(pov,"  direction < 0, 0, %.6g >\n",DIRECTION);
  fprintf(pov,"  look_at < 0, 0, %.6g >\n",eyeDist);
  fputs(      "  sky < 0, 1, 0 >\n",pov);
  fprintf(pov,"  up < 0, %.6g, 0 >\n",WINDOW_UP);
  fprintf(pov,"  right < %.6g, 0, 0 >\n",WINDOW_RIGHT);
  
/*
 * Tricky: rotate the camera into the right direction line of sight vector.
 * The look_at vector (as above, initially along +Z) and the sky vector
 * (initially along +Y) must be correctly rotated into their defined values.
 * Then translate the camera into place. All this is necessary because POV
 * seems to do some weird things to move the camera instead of rotating it
 * when a "sky" is defined.
 */
    
  lineOfSight.x=LOOK_AT_X-eye.x;
  lineOfSight.y=LOOK_AT_Y-eye.y;
  lineOfSight.z=LOOK_AT_Z-eye.z;
  XZproj=SQRT((lineOfSight.x*lineOfSight.x)+(lineOfSight.z*lineOfSight.z));
  if (XZproj==0.0) {
    rot2_X= (lineOfSight.y>0.0) ? -M_PI_2 : M_PI_2;
    rot3_Y=0.0;
  } else {
    rot2_X=-ASIN(lineOfSight.y/eyeDist);
    rot3_Y= (lineOfSight.x>0.0) ?
        ACOS(lineOfSight.z/XZproj) : -ACOS(lineOfSight.z/XZproj);
  }
  invSky.x=SKY_X;
  invSky.y=SKY_Y;
  invSky.z=SKY_Z;
  Apply_Rotation(0.0,1.0,0.0,-rot3_Y,&(invSky.x),&(invSky.y),&(invSky.z));
  Apply_Rotation(1.0,0.0,0.0,-rot2_X,&(invSky.x),&(invSky.y),&(invSky.z));
  rot1_Z= (invSky.x>0.0) ? -ACOS(invSky.y) : ACOS(invSky.y);
  sprintf(tmp,"  rotate < 0, 0, %.7g >\n",DEGREES(rot1_Z));
  if (strcmp(tmp,"  rotate < 0, 0, 0 >\n")!=0)
    fprintf(pov,"%s",tmp);
  sprintf(tmp,"  rotate < %.7g, 0, 0 >\n",DEGREES(rot2_X));
  if (strcmp(tmp,"  rotate < 0, 0, 0 >\n")!=0)
    fprintf(pov,"%s",tmp);
  sprintf(tmp,"  rotate < 0, %.7g, 0 >\n",DEGREES(rot3_Y));
  if (strcmp(tmp,"  rotate < 0, 0, 0 >\n")!=0)
    fprintf(pov,"%s",tmp);
  fputs(      "  translate Camera_Position\n",pov);

  fprintf(pov,"}\n\n");

  Light.x=CAMERA_X+cameraDist*(SKY_X)*0.65+cameraDist*(unitRight.x)*0.4;
  Light.y=CAMERA_Y+cameraDist*(SKY_Y)*0.65+cameraDist*(unitRight.y)*0.4;
  Light.z=CAMERA_Z+cameraDist*(SKY_Z)*0.65+cameraDist*(unitRight.z)*0.4;
  fprintf(pov,"light_source { < %.5g, %.5g, %.5g > color rgb < 1, 1, 1 > }\n",
              Light.x,Light.y,Light.z);

  if (whichEye==MONO)
    fputs("light_source { Camera_Position",pov);
  else
    fprintf(pov,"light_source { < %.5g, %.5g, %.5g >",
                CAMERA_X,CAMERA_Y,CAMERA_Z);
  fputs(" color rgb < 0.5, 0.5, 0.5 > }\n\n",pov);
  
  fprintf(pov,"background { color rgb < %.5g, %.5g, %.5g > }\n\n",
              Background.r, Background.g, Background.b);
  
  fprintf(pov,"#declare %s =\n",FINISH);
  fputs(FINISH_DECLARATION,pov);
  
  puts("Writing color and atom definitions...");
  Write_Color_and_Atom_Declarations();
  if (model!=CPK) {
    puts("Writing bond definitions...");
    Write_Bond_Declarations();
  }
  fprintf(pov,"\n#include \"%s\"\n",incName);
  fclose(pov);
  printf("Pov header %s created.\n\n",fileName);
  return;
}


/*------- adjust view - rotate, translate, zoom, print view params -------*/

void Adjust_View(void)
{
  char choice, tmp[80], quit=9, changed=FALSE;
  F_TYPE num, distance;
  Vector translate;
  
  do {
    puts("");
    puts(  "Options: 1) Rotate in plane");
    puts(  "         2) Rotate around vertical axis");
    puts(  "         3) Rotate around horizontal axis");
    puts(  "         4) Translate horizontal");
    puts(  "         5) Translate vertical");
    puts(  "         6) Change distance to camera");
    puts(  "         7) Change zoom factor");
    printf("         %u) Print view parameters\n",quit-1);
    printf("         %u) Return",quit);
    printf("\nWhich [%u] ? ",quit);
    choice= (Get_User_Input(tmp)[0]=='\0') ? quit : atoi(tmp);
    if (choice<1 || choice>quit) choice=quit;
    puts("");
    
    if (choice==1) {
      printf("Positive value rotates clockwise. How many degrees [30] ? ");
      num= (Get_User_Input(tmp)[0]=='\0') ? 30.0 : atof(tmp);
      if (num!=0.0) {
        Apply_Rotation(LOOK_AT_X-CAMERA_X,LOOK_AT_Y-CAMERA_Y,
                       LOOK_AT_Z-CAMERA_Z,RADIANS(num),
                       &SKY_X,&SKY_Y,&SKY_Z);
        Apply_Rotation(LOOK_AT_X-CAMERA_X,LOOK_AT_Y-CAMERA_Y,
                       LOOK_AT_Z-CAMERA_Z,RADIANS(num),
                       &(unitRight.x),&(unitRight.y),&(unitRight.z));
        changed=TRUE;
      }
    } else if (choice==2) {
      printf("Positive value rotates left-to-right. How many degrees [30] ? ");
      num= (Get_User_Input(tmp)[0]=='\0') ? 30.0 : atof(tmp);
      if (num!=0.0) {
        Apply_Rotation(SKY_X,SKY_Y,SKY_Z,
                       RADIANS(num),
                       &(unitRight.x),&(unitRight.y),&(unitRight.z));
        translate.x=LOOK_AT_X-CAMERA_X;
        translate.y=LOOK_AT_Y-CAMERA_Y;
        translate.z=LOOK_AT_Z-CAMERA_Z;
        Apply_Rotation(SKY_X,SKY_Y,SKY_Z,
                       RADIANS(num),
                       &(translate.x),&(translate.y),&(translate.z));
        CAMERA_X=LOOK_AT_X-translate.x;
        CAMERA_Y=LOOK_AT_Y-translate.y;
        CAMERA_Z=LOOK_AT_Z-translate.z;
        changed=TRUE;
      }
    } else if (choice==3) {
      printf("Positive value rotates top-to-bottom. How many degrees [30] ? ");
      num= (Get_User_Input(tmp)[0]=='\0') ? 30.0 : atof(tmp);
      if (num!=0.0) {
        Apply_Rotation(unitRight.x,unitRight.y,unitRight.z,
                       RADIANS(num),
                       &SKY_X,&SKY_Y,&SKY_Z);
        translate.x=LOOK_AT_X-CAMERA_X;
        translate.y=LOOK_AT_Y-CAMERA_Y;
        translate.z=LOOK_AT_Z-CAMERA_Z;
        Apply_Rotation(unitRight.x,unitRight.y,unitRight.z,
                       RADIANS(num),
                       &(translate.x),&(translate.y),&(translate.z));
        CAMERA_X=LOOK_AT_X-translate.x;
        CAMERA_Y=LOOK_AT_Y-translate.y;
        CAMERA_Z=LOOK_AT_Z-translate.z;
        changed=TRUE;
      }
    } else if (choice==4) {
      printf("Positive value moves right. How many units [10] ? ");
      num= (Get_User_Input(tmp)[0]=='\0') ? 10.0 : atof(tmp);
      if (num!=0.0) {
        CAMERA_X-=num*unitRight.x;
        CAMERA_Y-=num*unitRight.y;
        CAMERA_Z-=num*unitRight.z;
        LOOK_AT_X-=num*unitRight.x;
        LOOK_AT_Y-=num*unitRight.y;
        LOOK_AT_Z-=num*unitRight.z;
        changed=TRUE;
     }
    } else if (choice==5) {
      printf("Positive value moves up. How many units [10] ? ");
      num= (Get_User_Input(tmp)[0]=='\0') ? 10.0 : atof(tmp);
      if (num!=0.0) {
        CAMERA_X-=num*SKY_X;
        CAMERA_Y-=num*SKY_Y;
        CAMERA_Y-=num*SKY_Z;
        LOOK_AT_X-=num*SKY_X;
        LOOK_AT_Y-=num*SKY_Y;
        LOOK_AT_Z-=num*SKY_Z;
        changed=TRUE;
      }
    } else if (choice==6) {
      distance=Distance(CAMERA_X,CAMERA_Y,CAMERA_Z,
                        LOOK_AT_X,LOOK_AT_Y,LOOK_AT_Z);
      printf("Current distance is %.5g. Change to [%.5g] ? ",distance,distance);
      num= (Get_User_Input(tmp)[0]=='\0') ? distance : atof(tmp);
      if (num!=distance) {
        CAMERA_X=LOOK_AT_X+(num/distance)*(CAMERA_X-LOOK_AT_X);
        CAMERA_Y=LOOK_AT_Y+(num/distance)*(CAMERA_Y-LOOK_AT_Y);
        CAMERA_Z=LOOK_AT_Z+(num/distance)*(CAMERA_Z-LOOK_AT_Z);
        changed=TRUE;
      }
    } else if (choice==7) {
      printf("Current zoom is %.5g. Change to [%.5g] ? ",DIRECTION,DIRECTION);
      num= (Get_User_Input(tmp)[0]=='\0') ? DIRECTION : atof(tmp);
      if (num!=DIRECTION) {
        DIRECTION=num;
        changed=TRUE;
      }
    } else if (choice==(quit-1)) {
      printf("CameraX = %.5g\n",CAMERA_X);
      printf("CameraY = %.5g\n",CAMERA_Y);
      printf("CameraZ = %.5g\n",CAMERA_Z);
      printf("LookAtX = %.5g\n",LOOK_AT_X);
      printf("LookAtY = %.5g\n",LOOK_AT_Y);
      printf("LookAtZ = %.5g\n",LOOK_AT_Z);
      printf("SkyX = %.6g\n",SKY_X);
      printf("SkyY = %.6g\n",SKY_Y);
      printf("SkyZ = %.6g\n",SKY_Z);
      printf("WindowDistance = %.5g\n",DIRECTION);
      printf("WindowHeight   = %.5g\n",WINDOW_UP);
      printf("WindowWidth    = %.5g\n",WINDOW_RIGHT);
    }
  } while (choice!=quit);
  if (changed) {
    Write_POV_Header(MONO,povName);
    if (stereoPairCreated) {
      Write_POV_Header(LEFT_EYE,leftName);
      Write_POV_Header(RIGHT_EYE,rightName);
    }
  }
  return;
}


/*-------------------- the main program body --------------------------*/

void main(int argc, char *argv[])
{
  char *pdbName, *configName=NULL;
  char *tgaName=NULL, *leftTga=NULL, *rightTga=NULL;
  U_INT povWidth, povHeight, nBonds, central, bondedTo, i, j, quit;

#define TMPSZ 512

  char tmp[TMPSZ], con[6], choice=0, defaultChoice=0;
  PDBatomInfoType info;
  FILE *testTga;

  printf("\nWelcome to PovChem! This is version %s. For general\n",VERSION);
  puts("instructions see http://ludwig.scs.uiuc.edu/~paul/Manual.html\n");

/*---------------- parse command-line options --------------------*/

  for (i=1;i<argc;) {
    if (strcmp(argv[i],"-default")==0) {
      if (i<(argc-1)) {
        if (strcmp(argv[i+1],"-config")==0) { i=argc+1; break; }
        if (strcmp(argv[i+1],"-table")==0) { i=argc+1; break; }
        strcpy(tmp,argv[i+1]);
        isInteractive=FALSE;
        i+=2;
      } else { i=argc+1; break; }
    } else if (strcmp(argv[i],"-config")==0) {
      if (i<(argc-1)) {
        if (strcmp(argv[i+1],"-default")==0) { i=argc+1; break; }
        if (strcmp(argv[i+1],"-table")==0) { i=argc+1; break; }
        New_String(&configName,argv[i+1]);
        i+=2;
      } else { i=argc+1; break; }
    } else if (strcmp(argv[i],"-table")==0) {
      if (i<(argc-1)) {
        if (strcmp(argv[i+1],"-default")==0) { i=argc+1; break; }
        if (strcmp(argv[i+1],"-config")==0) { i=argc+1; break; }
        New_String(&PERIODIC_TABLE_FILE,argv[i+1]);
        i+=2;
      } else { i=argc+1; break; }
    } else { i=argc+1; break; }
  }

  if (i==(argc+1) || argc>7) {
    puts("\nUsage:");
    puts(  "povchem [-default <pdb filename>] [-config <config filename>]");
    puts(  "        [-table <periodic table filename>]\n");
    exit(EXIT_FAILURE);
  }


/*----------- Parse periodic table and config file -------------*/

  Assign_Defaults(&configName); /* must be read first in case periodic table
                                   is in non-default location.              */
  Parse_Periodic_Table();
  

/*------------ Get filenames and open files ---------------*/

  if (isInteractive) {
    printf("\nInput file [.pdb] : ");
    Get_User_Input(tmp);
  }
  if (strstr(tmp,".pdb")==0) strcat(tmp,".pdb");
  New_String(&pdbName,tmp);
  if ( (pdb=fopen(pdbName,"r"))==0 ) {
    printf("Can't read input file \"%s\".\n",pdbName);
    perror("Reason");
    exit(EXIT_FAILURE);
  }
  
  New_String(&incName,pdbName);
  *strrchr(incName,'.')='\0';
  strcat(incName,".inc");
  New_String(&povName,pdbName);
  *strrchr(povName,'.')='\0';
  strcat(povName,".pov");

  if (POVRAY) {
    New_String(&tgaName,pdbName);
    *strrchr(tgaName,'.')='\0';
    strcat(tgaName,".tga");
  }
  
  if ( (inc=fopen(incName,"w"))==0 ) {
    printf("Can't write include file \"%s\"\n",incName);
    perror("Reason");
    fclose(pdb);
    exit(EXIT_FAILURE);
  } else {
    fprintf(inc,"// This file was created by PovChem version %s\n\n",VERSION);
  }

  if ( (pov=fopen(povName,"w"))==0 ) {
    printf("Can't write pov header file \"%s\"\n",povName);
    perror("Reason");
    fclose(inc);
    fclose(pdb);
    exit(EXIT_FAILURE);
  }
  fclose(pov);


/*--------------- Decide which model to make ----------------*/

  if (isInteractive) {
    puts("\nWhat type of model? 1) CPK");
    puts  ("                    2) Ball and Stick");
    puts  ("                    3) Cylinder");
    printf("Choice [%u] ? ",DEFAULT_MODEL);
    model=( (Get_User_Input(tmp)[0]=='\0') ? DEFAULT_MODEL : atoi(tmp) );
    if (model<1 || model>3) model=DEFAULT_MODEL;
  } else
    model=DEFAULT_MODEL;

  switch (model) {


/*----------------- ball and stick -------------------------*/

    case BALLnSTICK: case CYLINDER:

                               /*  BALLnSTICK options  */
      if (model==BALLnSTICK) {
        if (isInteractive) {
          puts("\nSet Radius of the Balls: 1) Constant");
          puts(  "                         2) Proportional to vdW radius");
          printf("Choice [%u] ? ",DEFAULT_RADIUS_TYPE);
          ballSizeModel=( (Get_User_Input(tmp)[0]=='\0')
                               ? DEFAULT_RADIUS_TYPE : atoi(tmp) );
          if (ballSizeModel<1 || ballSizeModel>2)
            ballSizeModel=DEFAULT_RADIUS_TYPE;
        } else
          ballSizeModel=DEFAULT_RADIUS_TYPE;

        if (ballSizeModel==CONSTANT) {
          if (isInteractive) {
            printf("\nRadius of the Balls (in Angstroms) [%.2f] ? ",
                   DEFAULT_BALL_RADIUS);
            ballRadius=( (Get_User_Input(tmp)[0]=='\0')
                        ? DEFAULT_BALL_RADIUS : atof(tmp) );
            if (ballRadius<=0.0) ballRadius=DEFAULT_BALL_RADIUS;
          } else
            ballRadius=DEFAULT_BALL_RADIUS;
        } else {
          if (isInteractive) {
            printf("\nProportionality constant? [%.2f] ? ",
                   DEFAULT_PROPORTIONALITY);
            ballRadius=( (Get_User_Input(tmp)[0]=='\0')
                        ? DEFAULT_PROPORTIONALITY : atof(tmp) );
            if (ballRadius<=0.0) ballRadius=DEFAULT_PROPORTIONALITY;
          } else
            ballRadius=DEFAULT_PROPORTIONALITY;
        }

        if (isInteractive) {
          puts("\nColor the Balls by: 1) Element");
          puts  ("                    2) By PDB fields");
          puts  ("                    3) All the same");
          printf("Choice [%u] ? ",DEFAULT_COLORMODEL);
          sphereColorModel=( (Get_User_Input(tmp)[0]=='\0')
                        ? DEFAULT_COLORMODEL : atoi(tmp) );
          if (sphereColorModel<1 || sphereColorModel>3)
            sphereColorModel=DEFAULT_COLORMODEL;
        } else
          sphereColorModel=DEFAULT_COLORMODEL;

        if (sphereColorModel==ALL_THE_SAME) {
          allSphereColor.r=DEFAULT_BALL_RED;
          allSphereColor.g=DEFAULT_BALL_GREEN;
          allSphereColor.b=DEFAULT_BALL_BLUE;
          if (isInteractive) {
            printf("\nWhat color [%.5g, %.5g, %.5g] ? ",DEFAULT_BALL_RED,
                   DEFAULT_BALL_GREEN,DEFAULT_BALL_BLUE);
            Parse_Colors(Get_User_Input(tmp),&allSphereColor);
          }
        } else if (sphereColorModel==BY_PDB)
          Make_Tag_List(&firstSphereTag,&nSphereTags);
        
        if (isInteractive) {
          printf("\nRadius of the Sticks (for single bonds) [%.2f] ? ",
                 DEFAULT_STICK_RADIUS);
          cylinderRadius=( (Get_User_Input(tmp)[0]=='\0')
                          ? DEFAULT_STICK_RADIUS : atof(tmp) );
          if (cylinderRadius<=0.0) cylinderRadius=DEFAULT_STICK_RADIUS;
        } else
          cylinderRadius=DEFAULT_STICK_RADIUS;
          
        if (isInteractive) {
          puts("\nMake Radius of the Sticks: 1) All the same");
          puts(  "                           2) Thicker for higher order bonds");
          printf("Choice [%u] ? ",DEFAULT_THICKNESS);
          cylinderSizeModel=( (Get_User_Input(tmp)[0]=='\0')
                             ? DEFAULT_THICKNESS : atoi(tmp) );
          if (cylinderSizeModel<1 || cylinderSizeModel>2)
            cylinderSizeModel=DEFAULT_THICKNESS;
        } else
          cylinderSizeModel=DEFAULT_THICKNESS;
        
        if (isInteractive) {
          puts("\nColor the Sticks by: 1) Element");
          puts  ("                     2) By PDB fields");
          puts  ("                     3) All the same");
          printf("Choice [%u] ? ",DEFAULT_STICK_COLORMODEL);
          cylinderColorModel=( (Get_User_Input(tmp)[0]=='\0')
                              ? DEFAULT_STICK_COLORMODEL : atoi(tmp) );
          if (cylinderColorModel<1 || cylinderColorModel>3)
            cylinderColorModel=DEFAULT_STICK_COLORMODEL;
        } else
          cylinderColorModel=DEFAULT_STICK_COLORMODEL;
          
        if (cylinderColorModel==BY_PDB) {
          if (sphereColorModel==BY_PDB) {
            tmp[0]='\0';
            if (isInteractive) {
              printf("\nUse the same PDB fields as the Balls [%s] ? ",
                     (DEFAULT_STICKS_LIKE_BALLS ? "yes" : "no"));
              Get_User_Input(tmp);
            }
            if (tmp[0]=='\0')
              tmp[0]=(DEFAULT_STICKS_LIKE_BALLS ? 'y' : 'n');
            if (tolower(tmp[0])=='y') {
              firstCylinderTag=firstSphereTag;
              nCylinderTags=nSphereTags;
              stickTagType=sphereTagType;
            } else
              Make_Tag_List(&firstCylinderTag,&nCylinderTags);
          } else
            Make_Tag_List(&firstCylinderTag,&nCylinderTags);
        } else if (cylinderColorModel==ALL_THE_SAME) {
          allCylinderColor.r=DEFAULT_STICK_RED;
          allCylinderColor.g=DEFAULT_STICK_GREEN;
          allCylinderColor.b=DEFAULT_STICK_BLUE;
          if (isInteractive) {
            printf("\nWhat color [%.5g, %.5g, %.5g] ? ",DEFAULT_STICK_RED,
                   DEFAULT_STICK_GREEN,DEFAULT_STICK_BLUE);
            Parse_Colors(Get_User_Input(tmp),&allCylinderColor);
          }
        }

      } else {                /*  CYLINDER options  */

        if (isInteractive) {
          printf("\nRadius of the Cylinders (in Angstroms) [%.2f] ? ",
                 DEFAULT_CYLINDER_RADIUS);
          cylinderRadius=ballRadius=( (Get_User_Input(tmp)[0]=='\0')
                                     ? DEFAULT_CYLINDER_RADIUS : atof(tmp) );
          if (cylinderRadius<=0.0)
            ballRadius=cylinderRadius=DEFAULT_CYLINDER_RADIUS;
        } else
          ballRadius=cylinderRadius=DEFAULT_CYLINDER_RADIUS;
          
        if (isInteractive) {
          puts("\nColor the Cylinders by: 1) Element");
          puts  ("                        2) By PDB fields");
          puts  ("                        3) All the same");
          printf("Choice [%u] ? ",DEFAULT_COLORMODEL);
          sphereColorModel=cylinderColorModel=( (Get_User_Input(tmp)[0]=='\0')
                                            ? DEFAULT_COLORMODEL : atoi(tmp) );
          if (sphereColorModel<1 || sphereColorModel>3)
            sphereColorModel=cylinderColorModel=DEFAULT_COLORMODEL;
        } else
          sphereColorModel=cylinderColorModel=DEFAULT_COLORMODEL;
          
        if (sphereColorModel==ALL_THE_SAME) {
          allSphereColor.r=DEFAULT_STICK_RED;
          allSphereColor.g=DEFAULT_STICK_GREEN;
          allSphereColor.b=DEFAULT_STICK_BLUE;
          if (isInteractive) {
            printf("\nWhat color [%.5g, %.5g, %.5g] ? ",DEFAULT_STICK_RED,
                   DEFAULT_STICK_GREEN,DEFAULT_STICK_BLUE);
            Parse_Colors(Get_User_Input(tmp),&allSphereColor);
          }
          allCylinderColor.r=allSphereColor.r;
          allCylinderColor.g=allSphereColor.g;
          allCylinderColor.b=allSphereColor.b;
        } else if (sphereColorModel==BY_PDB) {
          Make_Tag_List(&firstSphereTag,&nSphereTags);
          firstCylinderTag=firstSphereTag;
          nCylinderTags=nSphereTags;
        }
        ballSizeModel=cylinderSizeModel=CONSTANT;

      }
      
                         /* options for both */
                         
      if (cylinderColorModel==BY_ELEMENT || cylinderColorModel==BY_PDB) {
        if (isInteractive) {
          puts("\nChange stick colors: 1) In the center");
          puts  ("                     2) At a point propotional to vdW radius");
          puts  ("                     3) Shade smoothly");
          printf("Choice [%u] ? ",DEFAULT_BOND_SPLIT);
          stickSplit=( (Get_User_Input(tmp)[0]=='\0')
                        ? DEFAULT_BOND_SPLIT : atoi(tmp) );
          if (stickSplit<1 || stickSplit>3)
            stickSplit=DEFAULT_BOND_SPLIT;
        } else
          stickSplit=DEFAULT_BOND_SPLIT;
      }
      puts("");

/*------------------- count atoms, check for CONECT --------------*/

      nAtoms=0;
      nBonds=0;
      for ( ; fgets(tmp,TMPSZ,pdb)!=0 && strncmp(tmp,"END   ",6)!=0 ; ) {
        if ( strncmp(tmp,"ATOM  ",6)==0 || strncmp(tmp,"HETATM",6)==0)
          nAtoms++;
        if ( nBonds==0 && strncmp(tmp,"CONECT",6)==0 ) nBonds=1;
      }
      rewind(pdb);
      printf("Found %u atoms...",nAtoms);
      if ( !nBonds ) {
        puts("\n... but no CONECT (connectivity) infomation found!");
        puts("This is required for anything but CPK models. Aborting.");
        BOMB;
      }


/* -------------- read atoms into memory, write spheres ---------------*/

      if ( (Atoms=(Atom *)malloc(sizeof(Atom)*nAtoms))==NULL ) {
        printf("\nNot enough memory for atoms: %ux%u bytes.\n",
               (int)sizeof(Atom),nAtoms);
        BOMB;
      } else
        memAlloced+=nAtoms*sizeof(Atom);
      if ( model==BALLnSTICK && cylinderColorModel==BY_PDB &&
           (cylinderEndTags=
             (char **)malloc(sizeof(char *)*nAtoms))==NULL ) {
        puts("\nNot enough memory for cylinder info.");
        BOMB;
      } else
        memAlloced+=nAtoms*sizeof(char *);
      for ( i=0 ; i<nAtoms &&
              fgets(tmp,TMPSZ,pdb)!=0 && strncmp(tmp,"END   ",6)!=0 ; ) {
        if (strncmp(tmp,"ATOM  ",6)==0 || strncmp(tmp,"HETATM",6)==0) {
          Parse_PDB_Line(tmp,&info,FIND_ALL);
          X(i)=info.x;
          Y(i)=info.y;
          Z(i)=info.z;
          ZNUM(i)=info.Z;
          SERNUM(i)=info.serial;
          switch (sphereColorModel) {
            case BY_ELEMENT:
              strcpy(tmp,PeriodicTable[ZNUM(i)].name);
              break;
            case ALL_THE_SAME:
              if (ballSizeModel==CONSTANT)
                strcpy(tmp,"All");
              else
                strcpy(tmp,PeriodicTable[ZNUM(i)].name);
              break;
            case BY_PDB:
              New_String(&(TAG(i)),
                         Make_Tag_From_PDB_Info(info,&sphereColorModel));
              strcpy(tmp,TAG(i));
              break;
          }
          Write_Atom(X(i),Y(i),Z(i),ZNUM(i),tmp);
          if (model==BALLnSTICK && cylinderColorModel==BY_PDB) {
            New_String(&(cylinderEndTags[i]),
                      Make_Tag_From_PDB_Info(info,&cylinderColorModel));
          }
          i++;
        }
      }
      rewind(pdb);
      fflush(inc);
      printf("\n ...loaded into memory.\nWrote %u spheres...",i);

/*----------- read connectivity list into memory, write cylinders ---------*/

      if ( (BondsTo=(U_INT *)
              malloc(nAtoms*MAX_BONDS_PER_ATOM*sizeof(U_INT)))==NULL ||
           (BondOrders=(unsigned char *)
              malloc(nAtoms*MAX_BONDS_PER_ATOM*sizeof(unsigned char)))==NULL ) {
        printf("\nNot enough memory for bonds.\n");
        BOMB;
      } else
        memAlloced+=(nAtoms*MAX_BONDS_PER_ATOM*
                     (sizeof(U_INT)+sizeof(unsigned char)) );
      for (i=0;i<nAtoms;i++) {
        for (j=0;j<MAX_BONDS_PER_ATOM;j++) {
          BoundTo(i,j)=nAtoms;
          BondOrder(i,j)=ZERO;
        }
      }
      for ( nBonds=0 ;
            fgets(tmp,TMPSZ,pdb)!=0 && strncmp(tmp,"END   ",6)!=0 ; ) {
        if ( strncmp(tmp,"CONECT",6)==0 ) {
          strncpy(con,tmp+6,5);
          con[5]='\0';
          central=atoi(con);
          i=0;
          do {
            con[0]=con[1]=con[2]=con[3]=con[4]='\0';
            strncpy(con,tmp+11+i*5,5);
            if (strpbrk(con,"0123456789")==NULL)
              i=37;
            else {
              bondedTo=atoi(con);
              if (central==bondedTo)
                printf("Can't bind atom %u to itself! Bond ignored.\n",central);
              else {
                Add_Bond(central,bondedTo,SINGLE);
                nBonds++;
                i++;
              }
            }
          } while (i<4);
        }
      }
      printf("\n ...and read %u bonds into memory.\n",nBonds);

      fclose(pdb);

      tmp[0]='\0';
      if (isInteractive) {
        printf("\nDo you wish to search for H-Bonds [%s] ? ",
               DEFAULT_FIND_H_BONDS ? "yes" : "no");
        Get_User_Input(tmp);
      }
      if (tmp[0]=='\0')
        tmp[0]=(DEFAULT_FIND_H_BONDS ? 'y' : 'n');
        
      if (tolower(tmp[0])=='y') {
        
        HBondColor.r=DEFAULT_H_BOND_RED;
        HBondColor.g=DEFAULT_H_BOND_GREEN;
        HBondColor.b=DEFAULT_H_BOND_BLUE;

        if (isInteractive) {
          printf("Use default parameters [yes] ? ");
          Get_User_Input(tmp);
        } else
          tmp[0]='y';
        
        if (tolower(tmp[0])!='y' && tmp[0]!='\0') {
          printf("H-Bond radius [%.2f] ? ",DEFAULT_H_BOND_RADIUS);
          HBondRadius=( (Get_User_Input(tmp)[0]=='\0')
            ? DEFAULT_H_BOND_RADIUS : atof(tmp));
          if (HBondRadius<=0.0) HBondRadius=DEFAULT_H_BOND_RADIUS;
          printf("Number of dashes in a typical H-bond [%u] ? ",
                 DEFAULT_H_BOND_DASHES);
          HBondDashes=( (Get_User_Input(tmp)[0]=='\0')
            ? DEFAULT_H_BOND_DASHES : atof(tmp));
          if (HBondDashes<=0) HBondDashes=DEFAULT_H_BOND_DASHES;
          printf("What color [%.5g, %.5g, %.5g] ? ",DEFAULT_H_BOND_RED,
                 DEFAULT_H_BOND_GREEN,DEFAULT_H_BOND_BLUE);
          Parse_Colors(Get_User_Input(tmp),&HBondColor);
          printf("Donors [%s] ? ",H_BOND_DONORS);
          if (Get_User_Input(tmp)[0]!='\0') {
            New_String(&H_BOND_DONORS,tmp);
            Make_H_Bonder_List(H_BOND_DONORS);
          }
          printf("Acceptors [%s] ? ",H_BOND_ACCEPTORS);
          if (Get_User_Input(tmp)[0]!='\0') {
            New_String(&H_BOND_ACCEPTORS,tmp);
            Make_H_Bonder_List(H_BOND_ACCEPTORS);
          }
          printf("Minimum H..A distance? [%.2f] ? ",DEFAULT_MINDIST);
          HBondMinDist=( (Get_User_Input(tmp)[0]=='\0')
                         ? DEFAULT_MINDIST : atof(tmp));
          if (HBondMinDist<=0) HBondMinDist=DEFAULT_MINDIST;
          printf("Maximum H..A distance? [%.2f] ? ",DEFAULT_MAXDIST);
          HBondMaxDist=( (Get_User_Input(tmp)[0]=='\0')
                         ? DEFAULT_MAXDIST : atof(tmp));
          if (HBondMaxDist<=0) HBondMaxDist=DEFAULT_MAXDIST;
          printf("Minimum DH..A angle? [%.2f] ? ",DEFAULT_MIN_DHA_ANGLE);
          HBondMinDHA=( (Get_User_Input(tmp)[0]=='\0')
                         ? DEFAULT_MIN_DHA_ANGLE : atof(tmp));
          if (HBondMinDHA<=0) HBondMinDHA=DEFAULT_MIN_DHA_ANGLE;
          printf("Minimum H..AR angle? [%.2f] ? ",DEFAULT_MIN_HAR_ANGLE);
          HBondMinHAR=( (Get_User_Input(tmp)[0]=='\0')
                         ? DEFAULT_MIN_HAR_ANGLE : atof(tmp));
          if (HBondMinHAR<=0) HBondMinHAR=DEFAULT_MIN_HAR_ANGLE;
        } else {
          HBondRadius=DEFAULT_H_BOND_RADIUS;
          HBondDashes=DEFAULT_H_BOND_DASHES;
          HBondMinDist=DEFAULT_MINDIST;
          HBondMaxDist=DEFAULT_MAXDIST;
          HBondMinDHA=DEFAULT_MIN_DHA_ANGLE;
          HBondMinHAR=DEFAULT_MIN_HAR_ANGLE;
        }

        Find_Hydrogen_Bonds();
      }
      puts("");
      
      Write_Bonds();
      fclose(inc);
      free(BondsTo);
      free(BondOrders);
      free(Atoms);
      
      break;


/*--------------- CPK ---------------------*/

    case CPK:
    
      if (isInteractive) {
        puts("\nColor the Atoms by: 1) Element");
        puts  ("                    2) By PDB fields");
        puts  ("                    3) All the same");
        printf("Choice [%u] ? ",DEFAULT_COLORMODEL);
        sphereColorModel=( (Get_User_Input(tmp)[0]=='\0')
                          ? DEFAULT_COLORMODEL : atoi(tmp) );
        if (sphereColorModel<1 || sphereColorModel>3)
          sphereColorModel=DEFAULT_COLORMODEL;
      } else
        sphereColorModel=DEFAULT_COLORMODEL;

      if (sphereColorModel==ALL_THE_SAME) {
        allSphereColor.r=DEFAULT_BALL_RED;
        allSphereColor.g=DEFAULT_BALL_GREEN;
        allSphereColor.b=DEFAULT_BALL_BLUE;
        if (isInteractive) {
          printf("\nWhat color [%.5g, %.5g, %.5g] ? ",DEFAULT_BALL_RED,
                  DEFAULT_BALL_GREEN,DEFAULT_BALL_BLUE);
          Parse_Colors(Get_User_Input(tmp),&allSphereColor);
        }
      } else if (sphereColorModel==BY_PDB) {
        Make_Tag_List(&firstSphereTag,&nSphereTags);
      }


      puts("");

      for (nAtoms=0;fgets(tmp,TMPSZ,pdb)!=0 && strncmp(tmp,"END   ",6)!=0;) {
        if ( strncmp(tmp,"ATOM  ",6)==0 || strncmp(tmp,"HETATM",6)==0) {
          Parse_PDB_Line(tmp,&info,FIND_ALL);
          switch (sphereColorModel) {
            case BY_ELEMENT: case ALL_THE_SAME:
              strcpy(tmp,PeriodicTable[info.Z].name);
              break;
            case BY_PDB:
              strcpy(tmp,Make_Tag_From_PDB_Info(info,&sphereColorModel));
              break;
          }
          Write_Atom(info.x,info.y,info.z,info.Z,tmp);
          nAtoms++;
        }
      }
      printf("\nSuccessfully wrote %u spheres.\n",nAtoms);
      fclose(pdb);
      fclose(inc);
      break;
  }
  
  printf("\n%s created.\n",incName);
  
/*---- Now the .inc file is written; we need to set up the actual
       .pov file with appropriate color definitions, view, etc... -----*/
  
  aveX=aveX/((F_TYPE)nAtoms);
  if (aveX<0.0001 && aveX>-0.0001) aveX=0.0; 
  aveY=aveY/((F_TYPE)nAtoms);
  if (aveY<0.0001 && aveY>-0.0001) aveY=0.0; 
  aveZ=aveZ/((F_TYPE)nAtoms);      
  if (aveZ<0.0001 && aveZ>-0.0001) aveZ=0.0; 

  Set_Up_View();
  if (sphereColorModel==BY_PDB)
    Assign_Colors_To_Tags(&sphereColorModel);
  if (model==BALLnSTICK && cylinderColorModel==BY_PDB &&
      !DEFAULT_STICKS_LIKE_BALLS)
    Assign_Colors_To_Tags(&cylinderColorModel);
  Write_POV_Header(MONO,povName);

  quit=3;
  if (POVRAY) {
    quit+=2;
    if (VIEWER) quit++;
  }
  defaultChoice=quit;

  do {
    if (isInteractive) {
      if (POVRAY) {
        printf("\nNow you can: 1) Make");
        if (VIEWER) printf(" and view");
        puts(" a preview picture");
        printf(  "             2) Make a full-quality");
        if (stereoPairCreated)
          printf(" (stereo)");
        puts(" picture");
      }
      if (VIEWER)
        puts  (  "             3) View your picture");
      
      if (!POVRAY && !VIEWER)
        printf("\nNow you can: ");
      else
        printf(  "             ");
      printf("%u) Make stereo pair\n",quit-2);
      printf(    "             %u) Adjust the View\n",quit-1);
      printf(    "             %u) Quit\n",quit);
      printf("Which [%u] ? ",defaultChoice);
      choice=( (Get_User_Input(tmp)[0]=='\0') ? defaultChoice : atoi(tmp) );
      if (choice<1 || choice>quit) choice=defaultChoice;

    } else {

      if (choice==0 && DEFAULT_MAKE_STEREO_PAIR)
        choice=quit-2;
      else if (RUN_POVRAY && (choice==0 || choice==quit-2))
        choice=2;
      else
        choice=quit;
    }
     
    if (POVRAY && choice==1) {

      switch (imageOrientation) {
        case LANDSCAPE: default:
          povWidth=240;
          break;
        case PORTRAIT:
          povWidth=180;
          break;
        case USER_DEFINED:
          if (WINDOW_RIGHT>WINDOW_UP)
            povWidth=240;
          else
            povWidth=180;
          break;
      }
      povHeight=Round((WINDOW_UP/WINDOW_RIGHT)*povWidth);
      sprintf(tmp,"%s %s -A +W%u +H%u +I%s +O%s",
                  POVRAY,POV_OPTIONS,povWidth,povHeight,
                  povName,tgaName);
      system(tmp);

      if (VIEWER) {
        if ((testTga=fopen(tgaName,"r"))!=NULL) {
          fclose(testTga);
          printf("\n%s present! Viewing...\n",tgaName);
          sprintf(tmp,"%s %s",VIEWER,tgaName);
          system(tmp);
        } else {
          sprintf(tmp,"\nCan't open your picture '%s' because",tgaName);
          perror(tmp);
        }
      }


    } else if (POVRAY && choice==2) {

      povWidth=POV_WIDTH;
      if (isInteractive) {
        printf("\nHow wide [%u] ? ",povWidth);
        povWidth=((Get_User_Input(tmp)[0]=='\0') ? povWidth : atoi(tmp));
      }
      povHeight=Round((WINDOW_UP/WINDOW_RIGHT)*povWidth);
      if (stereoPairCreated) {
        sprintf(tmp,"%s %s +W%u +H%u +I%s +O%s",
                    POVRAY,POV_OPTIONS,povWidth,povHeight,
                    leftName,leftTga);
        system(tmp);
        sprintf(tmp,"%s %s +W%u +H%u +I%s +O%s",
                    POVRAY,POV_OPTIONS,povWidth,povHeight,
                    rightName,rightTga);
        system(tmp);
      } else {
        sprintf(tmp,"%s %s +W%u +H%u +I%s +O%s",
                    POVRAY,POV_OPTIONS,povWidth,povHeight,
                    povName,tgaName);
        system(tmp);
      }

    } else if (POVRAY && VIEWER && choice==3) {

      if ((testTga=fopen(tgaName,"r"))!=NULL) {
        fclose(testTga);
        printf("\n%s present! Viewing...\n",tgaName);
        sprintf(tmp,"%s %s",VIEWER,tgaName);
        system(tmp);
      } else {
        sprintf(tmp,"\nCan't open your picture '%s' because",tgaName);
        perror(tmp);
        puts("Did you forget to make the picture?");
      }

    } else if (choice==(quit-1)) {
      Adjust_View();
      
    } else if (choice==(quit-2)) {

      if (leftName==NULL) {
        strcpy(tmp,pdbName);
        *strrchr(tmp,'.')='\0';
        strcat(tmp,"-left.pov");
        New_String(&leftName,tmp);
        strcpy(tmp,pdbName);
        *strrchr(tmp,'.')='\0';
        strcat(tmp,"-right.pov");
        New_String(&rightName,tmp);
        if (POVRAY) {
          strcpy(tmp,leftName);
          *strrchr(tmp,'.')='\0';
          strcat(tmp,".tga");
          New_String(&leftTga,tmp);
          strcpy(tmp,rightName);
          *strrchr(tmp,'.')='\0';
          strcat(tmp,".tga");
          New_String(&rightTga,tmp);
        }
      }
      Write_POV_Header(LEFT_EYE,leftName);
      Write_POV_Header(RIGHT_EYE,rightName);
      stereoPairCreated=TRUE;
    }
    
  } while (choice!=quit);
  
/*
  printf("\nAllocated a total of about %iKb memory. Bye!\n\n",
         (int)(memAlloced/1024.0+1));
*/
  puts("");
  exit(EXIT_SUCCESS);
}