diff --git a/tools/pymol_asphere/README b/tools/pymol_asphere/README new file mode 100755 index 0000000000..df49f8d197 --- /dev/null +++ b/tools/pymol_asphere/README @@ -0,0 +1,21 @@ +asphere_vis: Tool for triangulating aspherical particles + to convert LAMMPS output to PyMol input + +Building the tool: + +cd src +make # may need to edit Makefile for your system + +This will place the asphere_vis executable in the bin directory. + +Documentation is in the doc directory. Note that the doc file says +the ellipsoid axes are specified by an "Ellipsoid" section in the data +file. In the current version of LAMPS, this is now a "Shapes" +section. + +See instructions on how to run an example ellipsoid viz in the +examples dir. + +Mike Brown, Sandia National Labs +wmbrown at sandia.gov +June 2007 diff --git a/tools/pymol_asphere/doc/asphere_vis.1 b/tools/pymol_asphere/doc/asphere_vis.1 new file mode 100755 index 0000000000..350fa57cff --- /dev/null +++ b/tools/pymol_asphere/doc/asphere_vis.1 @@ -0,0 +1,209 @@ +.if !'\*(.T'ps' .if !'\*(.T'html' .tm warning: eqn should have been given a `-T\*(.T' option +.if '\*(.T'html' .if !'ps'ps' .tm warning: eqn should have been given a `-Tps' option +.if '\*(.T'html' .if !'ps'ps' .tm warning: (it is advisable to invoke groff via: groff -Thtml -e) +.lf 1 /usr/share/groff/1.18.1.1/tmac/eqnrc +.\" Startup file for eqn. +.EQ +.nr 0C \n(.C +.cp 0 +.ds 10 +.cp \n(0C +.lf 63 +.EN +.lf 1 asphere_vis.manpage +.TH asphere_vis 1 "June 11, 2007" "asphere_vis (Graphics Utilities) 0.1" "Graphics Utilities" +.SH NAME +\fBasphere_vis\fR - Tools for LAMMPS ellipsoid trajectory visualization in PyMol. +.PD 2 +.SH VERSION +.PD 1 +Version 0.1 +.PD 2 +.SH SYNOPSIS +.PD 1 +.TP +\fBasphere_vis\fR input_data_file input_dump_file output_py_file [\fB-b\fR] [\fB-c\fR \fIcolor_file\fR] [\fB-d\fR] [\fB-f\fR \fImax_frame\fR] [\fB-h\fR] [\fB-i\fR \fIstart_frame\fR \fIskip\fR \fIend_frame\fR] [\fB-n\fR \fInotice_level\fR] [\fB-r\fR \fIellip_res\fR] [\fB-s\fR] +.br +.PD 2 +.SH DESCRIPTION +.PD 1 +Tool for LAMMPS trajectory visualization in PyMol. The input is a LAMMPS 'data' file or a 'in' file with ellipsoid semi-axes specified using the ellipsoid command. The trajectory is input from a 'dump' file that must be generated using a custom style with the following arguments in order: +.PD 0 +.PP +.PD 1 + +.PD 0 +.TP +.PP +.PD 1 +\fItag type x y z quatw quati quatj quatk\fR +.PD 0 +.PP +.PD 1 + +.PD 2 +.SH PARAMETERS +.PD 1 +.TP +\fB-b\fR +When used with \fB-s\fR, the option will number the filenames based on the frame number. By default, they are numbered consequtively from zero. +.TP +\fB-c\fR \fIcolor_file\fR +.PD 0 +.TP +.PP +.PD 1 +Color the atom_types and set transparency based on the input file. The input file contains a space delimited set sequence of the color for an atom followed by the alpha. The color should be the string name and the alpha should be between 0 and 1. +.TP +\fB-d\fR +Use a LAMMPS input file rather than a data file for extracting atom type information. The input filename is specified as \fIinput_data_file\fR. +.TP +\fB-f\fR \fImax_frame\fR +.PD 0 +.TP +.PP +.PD 1 +Do not write more than \fImax_frame\fR frames to the output file. +.TP +\fB-h\fR +Print out the man page for help +.TP +\fB-i\fR \fIstart_frame\fR \fIskip\fR \fIend_frame\fR +.PD 0 +.TP +.PP +.PD 1 +Render the specified frame interval inclusive between \fIstart_frame\fR and \fIend_frame\fR. \fIskip\fR gives the number of frames to \fIskip\fR between each rendered frame. A value of 0 outputs every frame between \fIstart_frame\fR and \fIend_frame\fR. The first frame in the dump file is frame 0. +.TP +\fB-n\fR \fInotice_level\fR +.PD 0 +.TP +.PP +.PD 1 +Set the degree of program output. Use: +.PD 0 +.PP +.PD 1 + +.PD 0 +.PP +.PD 1 + \fB-n\fR 0 No output +.PD 0 +.PP +.PD 1 + \fB-n\fR 10 Normal program output +.PD 0 +.PP +.PD 1 + \fB-n\fR 20 Parameters useful for reproducing the results +.PD 0 +.PP +.PD 1 + \fB-n\fR 30 All output +.TP +\fB-r\fR \fIellip_res\fR +.PD 0 +.TP +.PP +.PD 1 +Resolution of ellipsoids in trajectory. The number of triangles per ellipsoid is equal to 2*(\fIellip_res\fR^2). Default is 10. +.TP +\fB-s\fR +Output the results into separate .py files. The filename and extension for the output files is taken from \fIoutput_py_file\fR. +.PD 2 +.SH AVAILABLE COLORS +.PD 1 + +.PD 0 +.PP +.PD 1 + black +.PD 0 +.PP +.PD 1 + blue +.PD 0 +.PP +.PD 1 + brown +.PD 0 +.PP +.PD 1 + cmyk_blue +.PD 0 +.PP +.PD 1 + cmyk_marine +.PD 0 +.PP +.PD 1 + deep +.PD 0 +.PP +.PD 1 + forest +.PD 0 +.PP +.PD 1 + green +.PD 0 +.PP +.PD 1 + grey +.PD 0 +.PP +.PD 1 + hotpink +.PD 0 +.PP +.PD 1 + magenta +.PD 0 +.PP +.PD 1 + marine +.PD 0 +.PP +.PD 1 + orange +.PD 0 +.PP +.PD 1 + purple +.PD 0 +.PP +.PD 1 + red +.PD 0 +.PP +.PD 1 + slate +.PD 0 +.PP +.PD 1 + teal +.PD 0 +.PP +.PD 1 + wheat +.PD 0 +.PP +.PD 1 + white +.PD 0 +.PP +.PD 1 + yellow +.PD 0 +.PP +.PD 1 + +.PD 2 +.SH KNOWN BUGS +.PD 1 +Comments are not allowed at any point between a section header and the end of the contents for a section in either the data file or the input file. +.PD 2 +.SH AUTHORS +.PD 1 +W. Michael Brown diff --git a/tools/pymol_asphere/doc/asphere_vis.manpage b/tools/pymol_asphere/doc/asphere_vis.manpage new file mode 100755 index 0000000000..40991503e7 --- /dev/null +++ b/tools/pymol_asphere/doc/asphere_vis.manpage @@ -0,0 +1,196 @@ +.TH asphere_vis 1 "June 11, 2007" "asphere_vis (Graphics Utilities) 0.1" "Graphics Utilities" +.SH NAME +\fBasphere_vis\fR - Tools for LAMMPS ellipsoid trajectory visualization in PyMol. +.PD 2 +.SH VERSION +.PD 1 +Version 0.1 +.PD 2 +.SH SYNOPSIS +.PD 1 +.TP +\fBasphere_vis\fR input_data_file input_dump_file output_py_file [\fB-b\fR] [\fB-c\fR \fIcolor_file\fR] [\fB-d\fR] [\fB-f\fR \fImax_frame\fR] [\fB-h\fR] [\fB-i\fR \fIstart_frame\fR \fIskip\fR \fIend_frame\fR] [\fB-n\fR \fInotice_level\fR] [\fB-r\fR \fIellip_res\fR] [\fB-s\fR] +.br +.PD 2 +.SH DESCRIPTION +.PD 1 +Tool for LAMMPS trajectory visualization in PyMol. The input is a LAMMPS 'data' file or a 'in' file with ellipsoid semi-axes specified using the ellipsoid command. The trajectory is input from a 'dump' file that must be generated using a custom style with the following arguments in order: +.PD 0 +.PP +.PD 1 + +.PD 0 +.TP +.PP +.PD 1 +\fItag type x y z quatw quati quatj quatk\fR +.PD 0 +.PP +.PD 1 + +.PD 2 +.SH PARAMETERS +.PD 1 +.TP +\fB-b\fR +When used with \fB-s\fR, the option will number the filenames based on the frame number. By default, they are numbered consequtively from zero. +.TP +\fB-c\fR \fIcolor_file\fR +.PD 0 +.TP +.PP +.PD 1 +Color the atom_types and set transparency based on the input file. The input file contains a space delimited set sequence of the color for an atom followed by the alpha. The color should be the string name and the alpha should be between 0 and 1. +.TP +\fB-d\fR +Use a LAMMPS input file rather than a data file for extracting atom type information. The input filename is specified as \fIinput_data_file\fR. +.TP +\fB-f\fR \fImax_frame\fR +.PD 0 +.TP +.PP +.PD 1 +Do not write more than \fImax_frame\fR frames to the output file. +.TP +\fB-h\fR +Print out the man page for help +.TP +\fB-i\fR \fIstart_frame\fR \fIskip\fR \fIend_frame\fR +.PD 0 +.TP +.PP +.PD 1 +Render the specified frame interval inclusive between \fIstart_frame\fR and \fIend_frame\fR. \fIskip\fR gives the number of frames to \fIskip\fR between each rendered frame. A value of 0 outputs every frame between \fIstart_frame\fR and \fIend_frame\fR. The first frame in the dump file is frame 0. +.TP +\fB-n\fR \fInotice_level\fR +.PD 0 +.TP +.PP +.PD 1 +Set the degree of program output. Use: +.PD 0 +.PP +.PD 1 + +.PD 0 +.PP +.PD 1 + \fB-n\fR 0 No output +.PD 0 +.PP +.PD 1 + \fB-n\fR 10 Normal program output +.PD 0 +.PP +.PD 1 + \fB-n\fR 20 Parameters useful for reproducing the results +.PD 0 +.PP +.PD 1 + \fB-n\fR 30 All output +.TP +\fB-r\fR \fIellip_res\fR +.PD 0 +.TP +.PP +.PD 1 +Resolution of ellipsoids in trajectory. The number of triangles per ellipsoid is equal to 2*(\fIellip_res\fR^2). Default is 10. +.TP +\fB-s\fR +Output the results into separate .py files. The filename and extension for the output files is taken from \fIoutput_py_file\fR. +.PD 2 +.SH AVAILABLE COLORS +.PD 1 + +.PD 0 +.PP +.PD 1 + black +.PD 0 +.PP +.PD 1 + blue +.PD 0 +.PP +.PD 1 + brown +.PD 0 +.PP +.PD 1 + cmyk_blue +.PD 0 +.PP +.PD 1 + cmyk_marine +.PD 0 +.PP +.PD 1 + deep +.PD 0 +.PP +.PD 1 + forest +.PD 0 +.PP +.PD 1 + green +.PD 0 +.PP +.PD 1 + grey +.PD 0 +.PP +.PD 1 + hotpink +.PD 0 +.PP +.PD 1 + magenta +.PD 0 +.PP +.PD 1 + marine +.PD 0 +.PP +.PD 1 + orange +.PD 0 +.PP +.PD 1 + purple +.PD 0 +.PP +.PD 1 + red +.PD 0 +.PP +.PD 1 + slate +.PD 0 +.PP +.PD 1 + teal +.PD 0 +.PP +.PD 1 + wheat +.PD 0 +.PP +.PD 1 + white +.PD 0 +.PP +.PD 1 + yellow +.PD 0 +.PP +.PD 1 + +.PD 2 +.SH KNOWN BUGS +.PD 1 +Comments are not allowed at any point between a section header and the end of the contents for a section in either the data file or the input file. +.PD 2 +.SH AUTHORS +.PD 1 +W. Michael Brown diff --git a/tools/pymol_asphere/doc/asphere_vis.pdf b/tools/pymol_asphere/doc/asphere_vis.pdf new file mode 100755 index 0000000000..df01a58669 Binary files /dev/null and b/tools/pymol_asphere/doc/asphere_vis.pdf differ diff --git a/tools/pymol_asphere/examples/README b/tools/pymol_asphere/examples/README new file mode 100644 index 0000000000..70ad3ee90d --- /dev/null +++ b/tools/pymol_asphere/examples/README @@ -0,0 +1,16 @@ +The asphere_vis tool needs an input script and a dump file +from a LAMMPS run. + +In examples/ellipse, you can run LAMMPS with in.ellipse +to generate dump.ellipse + +Copy in.ellipse and dump.ellipse to this dir. + +Run the tool to create a PyMol input file called ellipse.py by typing: + +../bin/asphere_vis in.ellipse dump.ellipse ellipse.py -c colors.ellipse -d + +Launch PyMol on this input and you should see a nice viz of +4 ellipsoids in a box of LJ particles: + +pymol ellipse.py diff --git a/tools/pymol_asphere/examples/colors.ellipse b/tools/pymol_asphere/examples/colors.ellipse new file mode 100644 index 0000000000..73fdefdd4e --- /dev/null +++ b/tools/pymol_asphere/examples/colors.ellipse @@ -0,0 +1,3 @@ +marine 0.3 +red 1.0 + diff --git a/tools/pymol_asphere/src/Makefile b/tools/pymol_asphere/src/Makefile new file mode 100755 index 0000000000..d7553f9238 --- /dev/null +++ b/tools/pymol_asphere/src/Makefile @@ -0,0 +1,200 @@ +#*************************************************************************** +# Makefile +# ------------------- +# +# _________________________________________________________________________ +# Build for the Graphics Utilities +# _________________________________________________________________________ +# +# begin : Thu June 9 2005 +# copyright : (C) 2003 by W. Michael Brown +# email : wmbrown@sandia.gov +# ***************************************************************************/ + +#Compiler type +#COMPILER = intel +COMPILER = gnu +#COMPILER = mpi +#COMPILER = mingw + +#Locations of outside objects relative to a source directory +HOBJ_DIR = ../obj + +BIN_DIR = ../bin + +ALL_DIR = . +ALL_LIB = $(HOBJ_DIR)/liball.a + +GRPHICS_DIR = . +GRPHICS_LIB = $(HOBJ_DIR)/libgraphics.a + +GRID_DIR = . +GRID_LIB = $(HOBJ_DIR)/libgrid.a + +MOL_DIR = . +MOL_LIB = $(HOBJ_DIR)/libmol.a + +MATH_DIR = . +MATH_LIB = $(HOBJ_DIR)/libmath.a + +EVERY_LIB = $(MOL_LIB) $(GRID_LIB) $(GRPHICS_LIB) $(MATH_LIB) $(ALL_LIB) + +MOLSIM_DIR = . +MOLSIM_LIB = $(HOBJ_DIR)/molsim.o $(HOBJ_DIR)/dynmif.o $(HOBJ_DIR)/dynmifq.o + +# Include directories +INC = -I$(ALL_DIR) -I$(MOLSIM_DIR) -I$(MOL_DIR) -I$(MATH_DIR) -I$(GRID_DIR) -I$(GRPHICS_DIR) + +ifeq ($(COMPILER),intel) + CPP = icpc # C++ Compiler + CC = icc # C compiler + AR = xiar #ar + DBUG = -g -DDEBUG -DNANCHECK #-Wall #-ansi + OPT = -O2 -xP -ipo -no-prec-div -static +endif + +ifeq ($(COMPILER),gnu) + CPP = g++ # C++ Compiler + CC = gcc # C compiler + AR = ar + DBUG = -g -DDEBUG -DNANCHECK -Wall -pedantic #-ansi + OPT = #-O3 +endif + +ifeq ($(COMPILER),mpi) + CPP = mpic++ -DMUSE_MPI # C++ Compiler + CC = mpicc -DMUSE_MPI # C compiler + AR = ar + DBUG = -g -DDEBUG -DNANCHECK #-Wall #-pedantic #-ansi + OPT = #-xN -O3 #-ipo -no-prec-div -static #-O3 +endif + +ifeq ($(COMPILER),mingw) + CPP = /cygdrive/c/MINGW/bin/g++ # C++ Compiler + CC = /cygdrive/c/MINGW/bin/gcc # C compiler + AR = /cygdrive/c/MINGW/bin/ar + DBUG = #-g -DDEBUG -DNANCHECK -Wall -pedantic #-ansi + OPT = -O3 -static +endif + +# Large file support? +LFSC = #-D_LARGEFILE_SOURCE `getconf LFS_CFLAGS` +LFSL = #`getconf LFS_LDFLAGS` `getconf LFS_LIBS` + +# GNU Scientific Library? +GSLC = #-DUSEGSL -I/usr/local/include/ +GSLL = #-lgsl -lgslcblas + +# GA Library? +LIBGAC = #-DLIBGA -I/usr/local/include/ +LIBGAL = #-lga -L/usr/local/lib/ -Wl,--allow-multiple-definition + +# Movie frame support? +MOVIE = #-DMOVIE + +# VTK ? +VTKH = #-DUSEVTK -I/usr/local/include/vtk-5.0 +VTKL = #-lvtkWidgets -lvtkHybrid -lvtkVolumeRendering -lvtkRendering -lvtkIO -lvtkGenericFiltering -lvtkGraphics -lvtkImaging -l vtkFiltering -lvtkCommon -L/usr/X11R6/lib/ -lGL -lXt -lSM -lICE -lX11 -lXext -lpthread -ldl + +CFLAGS = $(OPT) $(MOVIE) $(DBUG) $(INC) $(GSLC) $(LIBGAC) $(VTKH) -c +LFLAGS = $(OPT) +LLIBS = $(GSLL) $(LIBGAL) $(LFSL) $(VTKL) + +# Distribution Directories +DIST_BIN = /home/wmbrown/distbin/ +DIST_MAN = /home/wmbrown/cpp/manpages/man1/ +DIST_DOC = /home/wmbrown/cpp/doc/ + +OBJ_DIR = $(HOBJ_DIR) + +# Objects for this project +THIS_OBJ = $(OBJ_DIR)/asphere_vis.o $(GRPHICS_LIB) $(MATH_LIB) $(ALL_LIB) +EXECS = $(BIN_DIR)/asphere_vis + +all: $(EXECS) + +libraries: + cd $(ALL_DIR); make; cd $(MATH_DIR); make; cd $(GRID_DIR); make; \ + cd $(GRPHICS_DIR); make; + +ALL_OBJS = $(OBJ_DIR)/error.o $(OBJ_DIR)/commandline.o \ + $(OBJ_DIR)/misc.o + +$(OBJ_DIR)/error.o: error.h error.cpp + $(CPP) $(CFLAGS) -o $@ error.cpp + +$(OBJ_DIR)/commandline.o: commandline.h commandline.cpp + $(CPP) $(CFLAGS) -o $@ commandline.cpp + +$(OBJ_DIR)/misc.o: misc.h misc.cpp + $(CPP) $(CFLAGS) -o $@ misc.cpp + +$(ALL_LIB): $(ALL_OBJS) + $(AR) -crusv $(ALL_LIB) $(ALL_OBJS) + +GRPHICS_O = $(OBJ_DIR)/colors.o $(OBJ_DIR)/glsurface.o + +$(OBJ_DIR)/colors.o: colors.h colors.cpp + $(CPP) $(CFLAGS) -o $@ colors.cpp + +$(OBJ_DIR)/glsurface.o: glsurface.h glsurface.cpp + $(CPP) $(CFLAGS) -o $@ glsurface.cpp + +$(GRPHICS_LIB): $(GRPHICS_O) + $(AR) -crusv $(GRPHICS_LIB) $(GRPHICS_O) + +MATH_OBJS = $(OBJ_DIR)/cartesian.o $(OBJ_DIR)/miscm.o \ + $(OBJ_DIR)/spherical.o + +$(OBJ_DIR)/miscm.o: miscm.h miscm.cpp + $(CPP) $(CFLAGS) -o $@ miscm.cpp + +$(OBJ_DIR)/cartesian.o: cartesian.h cartesian.cpp + $(CPP) $(CFLAGS) -o $@ cartesian.cpp + +$(OBJ_DIR)/spherical.o: spherical.h spherical.cpp + $(CPP) $(CFLAGS) -o $@ spherical.cpp + +$(MATH_LIB): $(MATH_OBJS) + $(AR) -crusv $(MATH_LIB) $(MATH_OBJS) + +$(OBJ_DIR)/asphere_vis.o: asphere_vis.cpp + $(CPP) $(CFLAGS) -o $@ asphere_vis.cpp + +$(BIN_DIR)/asphere_vis: $(THIS_OBJ) + $(CPP) $(LFLAGS) -o $@ $(THIS_OBJ) $(LLIBS) +# +# Documentation +# +manpages: all + /bin/tcsh make_manpages.sh + +# +# Create a .tar distribution file +# +dist: all manpages + /bin/tcsh makedistribution.sh + + +# +# INSTALL to Mike's Directories +install: all manpages + /bin/cp $(EXECS) $(DIST_BIN); \ + /bin/cp ./manpages/*.1 $(DIST_MAN); \ + /bin/cp ./manpages/*.pdf $(DIST_DOC) + + +# +# Remove objects, cores, etc. +# + +clean: + rm -rf $(EXECS) $(THIS_OBJ) $(ALL_OBJ) $(MATH_OBJ) $(GRPHICS_O) + cd $(OBJ_DIR); rm -f *.o + +veryclean: clean + rm -rf *~ ./api ./manpages + +cleanproject: clean + cd $(ALL_DIR); make clean; cd $(MATH_DIR); make clean; \ + cd $(GRPHICS_DIR); make clean; diff --git a/tools/pymol_asphere/src/asphere_vis.cpp b/tools/pymol_asphere/src/asphere_vis.cpp new file mode 100644 index 0000000000..69bd353672 --- /dev/null +++ b/tools/pymol_asphere/src/asphere_vis.cpp @@ -0,0 +1,438 @@ +/*************************************************************************** + asphere_vis.cpp + ------------------- + + Convert a Lammps ellipsoid trajectory to a Pymol CGO trajectory + + __________________________________________________________________________ + This file is part of the Graphics Utilities package for command-line + access to Graphics Library functions +__________________________________________________________________________ + + begin : Fri Jan 12 2007 + copyright : (C) 2007 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#include "commandline.h" +#include "glsurface.h" +#include + +// Describe the program parameters +void Describe(CommandLine &cl,ostream &out); +// Parse the command line parameters +void HandleArgs(CommandLine &cl, int argc, char *argv[], Error *error); +// Parse past an ITEM line in the lammps dump file +bool parse_to(const char *token,ifstream &in); + +int main(int argc, char *argv[]) { + CommandLine cl; + Error error; + Colors colors; + FileIterator fi; + + // Parse the command line + HandleArgs(cl,argc,argv,&error); + + // ----------------- Set up file names + if (cl['s']) { + fi.set_file_header(a::namewoext(cl.argstring(' ',2))+"."); + fi.set_file_extensions("."+a::extension(cl.argstring(' ',2))); + fi.set_lead_zeros(4); + } + + // ----------------- Get the frame interval + unsigned start_frame=0; + unsigned skip=0; + unsigned end_frame=std::numeric_limits::max(); + if (cl['i']) { + start_frame=cl.argint('i',0); + skip=cl.argint('i',1); + end_frame=cl.argint('i',2); + } + + // ----------------- Get the atom type info from a data file + unsigned atom_types=0; + vector shapes; + ifstream in; + if (!cl['d']) { + a::fileopen(in,cl.argstring(' ',0),error); + while (!in.eof()) { + char iline[500]; + in.getline(iline,500); + vector tokens; + a::get_tokens(iline,tokens); + if (tokens.size()>2) + if (tokens[1]=="atom" && tokens[2]=="types") { + atom_types=atoi(tokens[0].c_str()); + break; + } + if (!in) + break; + } + if (atom_types==0) + error.generate_error(0,"asphere_vis", + "Could not find number of atom types in data file."); + // ----------------- Get the atom type shapes + if (parse_to("Shapes",in)) + for (unsigned i=0; i> tempd >> shape; + shape.x *= 0.5; + shape.y *= 0.5; + shape.z *= 0.5; + if (in.eof() || !in) + break; + shapes.push_back(shape); + char iline[500]; + in.getline(iline,500); + } + if (shapes.size()!=atom_types) { + error.buffer() << "Error reading shapes from Pair Coeffs section of " + << "data file. Read " << shapes.size() << " valid shapes, " + << "but expected " << atom_types; + error.addbuf(0,"asphere_vis"); + } + a::fileclose(in,error); + } else { + // ----------------- Get the atom type info from a input file + a::fileopen(in,cl.argstring(' ',0),error); + while (!in.eof()) { + string token; + in >> token; + if (token=="create_box") { + in >> atom_types; + shapes.assign(atom_types,cPt(0.5,0.5,0.5)); + } else if (token=="shape") { + unsigned type; + cPt shape; + in >> type >> shape; + shape.x *= 0.5; + shape.y *= 0.5; + shape.z *= 0.5; + if (type>atom_types) { + error.buffer() << "Error reading shapes from LAMMPS input file. " + << "I thought there were " << atom_types + << " atom types. But found an shape command for " + << "atom type: " << type; + error.addbuf(0,"asphere_vis"); + } + shapes[type-1]=shape; + } else { + char iline[500]; + in.getline(iline,500); + } + + if (!in && !in.eof()) + error.generate_error(0,"asphere_vis", + "Error reading from LAMMPS input file"); + } + a::fileclose(in,error); + } + if (atom_types==0) + error.generate_error(0,"asphere_vis","Found 0 atom types!"); + + // ----------------- Get the colors and alpha values for atom types + vector color_list; + vector alpha_list; + if (cl['c']) { + a::fileopen(in,cl.argstring('c',0),error); + for (unsigned i=0; i> color >> alpha; + if (in.eof() || !in) + error.generate_error(0,"asphere_vis", + "Improperly formatted color file."); + color_list.push_back(colors[color]); + alpha_list.push_back(alpha); + } + a::fileclose(in,error); + } + + a::fileopen(in,cl.argstring(' ',1),error); + ofstream out; + if (!cl['s']) { + a::fileopen(out,cl.argstring(' ',2),error); + } + + + // ----------------- Get the atom count + unsigned atom_count; + if (parse_to("ITEM: NUMBER OF ATOMS",in)) + in >> atom_count; + else + error.generate_error(0,"asphere_vis", + "Could not find ITEM: NUMBER OF ATOMS in input file."); + if (!in) + error.generate_error(0,"asphere_vis", + "Error reading ITEM: NUMBER OF ATOMS in input file."); + + // ----------------- Get the triangles per ellipsoid + unsigned ellip_res=10; + if (cl['r']) + ellip_res=cl.argint('r',0); + if (ellip_res==0) { + error.addwarning(0,9,"asphere_vis","Cannot use -r 0. Setting to 10."); + ellip_res=10; + } + + // ----------------- Get the bounding box + bool bound_found=false; + if (!cl['s']) { + if (parse_to("ITEM: BOX BOUNDS",in)) { + bound_found=true; + cPt bound[2]; + in >> bound[0].x >> bound[1].x; + in >> bound[0].y >> bound[1].y; + in >> bound[0].z >> bound[1].z; + GLSurface gls; + Vertex v; + v.transparency=1; + v.valid_normal=false; + v.color=colors["white"]; + for (unsigned i=0; i<2; i++) + for (unsigned j=0; j<2; j++) + for (unsigned k=0; k<2; k++) { + v.cpt=cPt(bound[i].x,bound[j].y,bound[k].z); + gls.addvertex(v); + } + gls.addline(0,1); + gls.addline(0,2); + gls.addline(0,4); + gls.addline(1,3); + gls.addline(1,5); + gls.addline(2,3); + gls.addline(2,6); + gls.addline(3,7); + gls.addline(4,5); + gls.addline(4,6); + gls.addline(5,7); + gls.addline(6,7); + gls.writelines(out,"gridb"); + out << "cmd.set('cgo_dot_width',8)\n"; + } else + error.addwarning(0,9,"asphere_vis", + "Could not find ITEM: BOX BOUNDS in input file. No box output."); + } + + if (!in) + error.generate_error(0,"asphere_vis", + "Error reading ITEM: BOX BOUNDS."); + + // ----------------- Generate the frames + unsigned frame=0; + unsigned max_frame=std::numeric_limits::max(); + if (cl['f']) + max_frame=cl.argint('f',0); + + colorPt color=colors["marine"]; + double alpha=1.0; + Vertex v; + v.color=color; + v.transparency=1.0; + + // ----------------- Get to the start frame + while (frameend_frame) + break; + if (!parse_to("ITEM: ATOMS",in)) + break; + GLSurface gls; + for (unsigned i=0; i> id >> atom_type; + cPt atom_center; + in >> atom_center; + Quaternion q; + in >> q; + if (!in) { + error.addwarning(0,9,"asphere_vis","Error reading frame: "+ + a::itoa(frame)); + break; + } + cPt radius(shapes[atom_type-1]); + if (radius.x == radius.y && radius.y == radius.z) { + v.cpt=atom_center; + if (cl['c']) { + v.color=color_list[atom_type-1]; + v.transparency=alpha_list[atom_type-1]; + } + gls.addvertex(v); + gls.add_sphere(gls.size_vertices()-1,radius.x); + } else { + if (cl['c']) { + color=color_list[atom_type-1]; + alpha=alpha_list[atom_type-1]; + } + gls.add_ellipsoid(atom_center,radius,q, + color,alpha,ellip_res); + } + } + if (!in) + break; + if (cl['s']) { + if (cl['b']) + fi.set_file_num(frame); + a::fileopen(out,fi.nextfilename(),error); + } + gls.writetris(out,"ellipse"); + if (gls.size_spheres()!=0) + gls.writespheres(out,"spheres"); + if (cl['s']) + a::fileclose(out,error); + wrote++; + frame++; + if (frame==max_frame) + break; + for (unsigned i=0; i " << name + << ".1' to generate a man page for this\n" + << "program and type 'man ./" << name << ".1' for help\n" + << "______________________________________________________________________\n"; + return; +} + +void HandleArgs(CommandLine &cl, int argc, char *argv[], Error *error) { + // Arguments + cl.addmanditory(' ',3); + cl.addargname(' ',"input_data_file"); + cl.addargname(' ',"input_dump_file"); + cl.addargname(' ',"output_py_file"); + cl.add('d',0); + cl.adddescription('d',"Use a LAMMPS input file rather than a data file for extracting atom type information. The input filename is specified as input_data_file."); + cl.add('f',1); + cl.addargname('f',"max_frame"); + cl.adddescription('f',"Do not write more than max_frame frames to the output file."); + cl.add('r',1); + cl.addargname('r',"ellip_res"); + cl.adddescription('r',"Resolution of ellipsoids in trajectory. The number of triangles per ellipsoid is equal to 2*(ellip_res^2). Default is 10."); + cl.add('c',1); + cl.addargname('c',"color_file"); + cl.adddescription('c',"Color the atom_types and set transparency based on the input file. The input file contains a space delimited set sequence of the color for an atom followed by the alpha. The color should be the string name and the alpha should be between 0 and 1."); + cl.add('s',0); + cl.adddescription('s',"Output the results into separate .py files. The filename and extension for the output files is taken from output_py_file."); + cl.add('i',3); + cl.addargname('i',"start_frame"); + cl.addargname('i',"skip"); + cl.addargname('i',"end_frame"); + cl.adddescription('i',"Render the specified frame interval inclusive between start_frame and end_frame. skip gives the number of frames to skip between each rendered frame. A value of 0 outputs every frame between start_frame and end_frame. The first frame in the dump file is frame 0."); + cl.add('b',0); + cl.adddescription('b',"When used with -s, the option will number the filenames based on the frame number. By default, they are numbered consequtively from zero."); + + // Stuff for every executable + cl.addhelp('h',0); + cl.adddescription('h',"Print out the man page for help"); + cl.add('n',1); + cl.addargname('n',"notice_level"); + cl.adddescription('n',"Set the degree of program output. Use: \n\n\t-n 0\tNo output\n\t-n 10\tNormal program output\n\t-n 20\tParameters useful for reproducing the results\n\t-n 30\tAll output"); + + // Short Description + cl.addtoman_chapter("NAME","Tools for LAMMPS ellipsoid trajectory visualization in PyMol."); + + // Version + cl.addtoman_chapter("VERSION","Version 0.1"); + + // Full Description + const string desc[5]={ + "Tool for LAMMPS trajectory visualization in PyMol. The input is a LAMMPS ", + "'data' file or a 'in' file with ellipsoid semi-axes specified using the ", + "ellipsoid command. The trajectory is input from a 'dump' file that must ", + "be generated using a custom style with the following arguments in order:\n", + ".TP\\fItag type x y z quatw quati quatj quatk\\fR\n" + }; + cl.addtoman_chapter("DESCRIPTION",5,desc); + + Colors colors; + cl.addtoman_chapter("AVAILABLE COLORS",colors.colorlist()); + + // bugs + const string bugs[3]={ + "Comments are not allowed at any point between a section header and ", + "the end of the contents for a section in either the data file or ", + "the input file." + }; + cl.addtoman_chapter("KNOWN BUGS",3,bugs); + + + // Authors + cl.addtoman_chapter("AUTHORS","W. Michael Brown"); + + // Parse the commandline + if (!cl.parse(argc,argv,error)) { + Describe(cl,cout); + error->generate_error(0,a::filenameonly(argv[0]),"Bad Command Line\n"); + } + + // Set the notice level + if (cl['n']) + error->note.set_notice_level(cl.argint('n',0)); + + // Generate a notice with the command line for records purposes + string cm=cl.program_name(); + for (int j=1; jnote.notice(19,"CommandLine",cm); + + // Output the help + if (cl['h']) { + cl.write_man_page(cout,"0.1","Graphics Utilities"); + exit(0); + } +} + +// Parse past an ITEM line in the lammps dump file +bool parse_to(const char * token,ifstream &in) { + char iline[5000]; + while (true) { + in.getline(iline,5000); + if (in.eof() || !in) + return false; + if (strcmp(token,iline)==0) + return true; + } +} diff --git a/tools/pymol_asphere/src/cartesian.cpp b/tools/pymol_asphere/src/cartesian.cpp new file mode 100755 index 0000000000..d412d2ee43 --- /dev/null +++ b/tools/pymol_asphere/src/cartesian.cpp @@ -0,0 +1,1093 @@ +/**************************************************************************** + * cartesian.cpp + * Shared stuff for dealing with Cartesian coordinates + * Also contains quaternion structures and operations + * + * Cartesian point of doubles: cPt + * Cartesian point of intergers: iPt + * Vector of doubles: vectorPt + * Color of doubles: colorPt + * Quaternion: Quaternion + * + * Can be accessed as cPt.x, cPt[X], colorPt[GREEN], iPt[I], etc. + * + * + * W. Michael Brown + * 5/26/03 + ****************************************************************************/ + +#include "cartesian.h" + +template class TwoD; +template class TwoD; +template class TwoD; +template ostream & operator<< (ostream &out, const TwoD &t); +template ostream & operator<< (ostream &out, const TwoD &t); +template istream & operator>> (istream &in, TwoD &t); +template istream & operator>> (istream &in, TwoD &t); +template class ThreeD; +template class ThreeD; +template class ThreeD; +template class ThreeD; +template ostream & operator<< (ostream &out, const ThreeD &t); +template istream & operator>> (istream &in, ThreeD &t); +template ostream & operator<< (ostream &out, const ThreeD &t); +template istream & operator>> (istream &in, ThreeD &t); +template ostream & operator<< (ostream &out, const ThreeD &t); +template istream & operator>> (istream &in, ThreeD &t); +template ThreeD operator+ (const unsigned one, + const ThreeD &two); +template ThreeD operator+(const double one, const ThreeD &two); +template ThreeD operator- (const unsigned one, + const ThreeD &two); +template ThreeD operator-(const double one, const ThreeD &two); +template ThreeD operator* (const unsigned one, + const ThreeD &two); +template ThreeD operator*(const double one, const ThreeD &two); +template ThreeD operator/ (const unsigned one, + const ThreeD &two); +template ThreeD operator/(const double one, const ThreeD &two); + + +// ------------------------ TwoD Stuff + +template +TwoD::TwoD() { +} + +// Assign both x and y the value +template +TwoD::TwoD(numtyp in) { + x=in; + y=in; +} + +// Assignment Constructor +template +TwoD::TwoD(numtyp cx, numtyp cy) { + x=cx; + y=cy; +} + +// Type conversion +template +TwoD::TwoD(const TwoD &two) { + x=numtyp(two.x); + y=numtyp(two.y); +} + +// Ball projection (onto xy-plane) +template +TwoD::TwoD(const Ball &ball) { + x=numtyp(cos(ball.theta)); + y=numtyp(sin(ball.theta)); +} + +template +numtyp & TwoD::operator[](unsigned i) { + switch(i) { + case X: return x; + case Y: return y; + } + return x; +} + +template +TwoD TwoD::operator + (const TwoD &two) const { + return TwoD(x+two.x, y+two.y); +} + +template +void TwoD::operator += (const TwoD &two) { + x+=two.x; + y+=two.y; +} + +template +TwoD TwoD::operator + (const numtyp two) const { + return TwoD(x+two, y+two); +} + +template +TwoD TwoD::operator - (const TwoD &two) const { + return TwoD(x-two.x, y-two.y); +} + +template +TwoD TwoD::operator - (const numtyp two) const { + return TwoD(x-two, y-two); +} + +template +TwoD TwoD::operator * (const numtyp two) const { + return TwoD(x*two, y*two); +} + +template +TwoD TwoD::operator * (const TwoD &two) const { + return TwoD(x*two.x,y*two.y); +} + +template +void TwoD::operator /= (const numtyp two) { + x/=two; + y/=two; +} + +template +bool TwoD::operator != (const TwoD &two) const { + if (x!=two.x || y!=two.y) + return true; + return false; +} + +// Move coordinates into array +template +void TwoD::to_array(numtyp *array) { + *array=x; + array++; + *array=y; +} + +// Set coordinates from array +template +void TwoD::from_array(numtyp *array) { + x=*array; + array++; + y=*array; +} + +// Dot Product +template +numtyp TwoD::dot(const TwoD &two) const { + return x*two.x+y*two.y; +} + +// Returns one of two normals to a line represented by vector +template +TwoD TwoD::normal() { + return TwoD(y,x*numtyp(-1)); +} + +template +numtyp TwoD::dist(const TwoD &two) const { + TwoD diff=*this-two; + return numtyp(sqrt(double(diff.dot(diff)))); +} + +template +numtyp TwoD::dist2(const TwoD &two) const { + TwoD diff=*this-two; + return diff.dot(diff); +} + +// Returns two +template +unsigned TwoD::dimensionality() { + return 2; +} + +// Returns the index of *this (for unsigned) in a square 3D array +/* s_size[0]=1Dsize */ +template +numtyp TwoD::array_index(vector &s_size) { + return y+x*s_size[0]; +} + +// Increment a 3D index from min to max (same as nested for) +template <> +bool TwoD::increment_index(TwoD &minp,TwoD &maxp) { + x++; + if (x!=maxp.x) + return true; + x=minp.x; + y++; + if (y!=maxp.y) + return true; + return false; +} + +// Return false if x or y is not within the inclusive range +template +bool TwoD::check_bounds(numtyp minp,numtyp maxp) { + if (xmaxp || y>maxp) + return false; + return true; +} + +template +ostream & operator<<(ostream &out, const TwoD &t) { + out << t.x << " " << t.y; + return out; +} + +template +istream & operator>>(istream &in, TwoD &t) { + in >> t.x >> t.y; + return in; +} + +// ------------------------ ThreeD Stuff + +// Assignment Constructor +template +ThreeD::ThreeD(numtyp cx, numtyp cy, numtyp cz) { + x=cx; y=cy; z=cz; +} + +// Assignment Constructor +template +ThreeD::ThreeD(numtyp in) { + x=in; y=in; z=in; +} + +// Type Conversion +template +ThreeD::ThreeD(const ThreeD& in) { + x=numtyp(in.x); y=numtyp(in.y); z=numtyp(in.z); +} + +// Type Conversion +template +ThreeD::ThreeD(const ThreeD& in) { + x=numtyp(in.x); y=numtyp(in.y); z=numtyp(in.z); +} + +// Type Conversion +template +ThreeD::ThreeD(const ThreeD& in) { + x=numtyp(in.x); y=numtyp(in.y); z=numtyp(in.z); +} + +// Type Conversion +template +ThreeD::ThreeD(const ThreeD& in) { + x=numtyp(in.x); y=numtyp(in.y); z=numtyp(in.z); +} + +// TwoD with (z-coordinate set to zero) +template +ThreeD::ThreeD(const TwoD& in) { + x=numtyp(in.x); y=numtyp(in.y); z=0; +} + +// Spherical Conversion +template +ThreeD::ThreeD(const Ball& ball) { + double sp=sin(ball.phi); + x=numtyp(sp*cos(ball.theta)); + y=numtyp(sp*sin(ball.theta)); + z=numtyp(cos(ball.phi)); +} + +// Spherical Conversion +template +ThreeD::ThreeD(const Ball& ball) { + double sp=sin(ball.phi); + x=numtyp(sp*cos(ball.theta)); + y=numtyp(sp*sin(ball.theta)); + z=numtyp(cos(ball.phi)); +} + +template +ThreeD::ThreeD() { +} + +template +numtyp & ThreeD::operator[](unsigned i) { + switch(i) { + case X: return x; + case Y: return y; + case Z: return z; + } + return x; +} + +template +numtyp ThreeD::operator[](unsigned i) const { + switch(i) { + case X: return x; + case Y: return y; + case Z: return z; + } + return x; +} + +template +void ThreeD::operator = (const ThreeD &two) { + #ifdef NANCHECK + //assert(a::not_nan(two.x) && a::not_nan(two.y) && a::not_nan(two.z)); + #endif + x=two.x; y=two.y; z=two.z; +} + +template +bool ThreeD::operator == (const ThreeD &two) const { + if (x!=two.x || y!=two.y || z!=two.z) + return false; + return true; +} + +template +bool ThreeD::operator != (const ThreeD &two) const { + if (x!=two.x || y!=two.y || z!=two.z) + return true; + return false; +} + +template +ThreeD ThreeD::operator + (const ThreeD &two) const { + return (ThreeD(x+two.x,y+two.y,z+two.z)); +} + +template +ThreeD ThreeD::operator + (const numtyp &two) const { + return (ThreeD(x+two,y+two,z+two)); +} + +template +void ThreeD::operator += (const numtyp &two) { + x+=two; y+=two; z+=two; +} + +template +void ThreeD::operator += (const ThreeD &two) { + x+=two.x; y+=two.y; z+=two.z; +} + +template +void ThreeD::operator -= (const ThreeD &two) { + x-=two.x; y-=two.y; z-=two.z; +} + +template +ThreeD ThreeD::operator - (const ThreeD &two) const { + return (ThreeD(x-two.x,y-two.y,z-two.z)); +} + +template +ThreeD ThreeD::operator - (const numtyp &two) const { + return (ThreeD(x-two,y-two,z-two)); +} + +template +void ThreeD::operator -= (const numtyp &two) { + x-=two; y-=two; z-=two; +} + +template +ThreeD ThreeD::operator * (const numtyp &two) const { + return (ThreeD(x*two,y*two,z*two)); +} + +template +ThreeD ThreeD::operator * (const ThreeD &two) const { + return ThreeD(x*two.x,y*two.y,z*two.z); +} + +template +ThreeD ThreeD::operator / (const ThreeD &two) const { + return ThreeD(x/two.x,y/two.y,z/two.z); +} + +// Dot Product +template +numtyp ThreeD::dot(const ThreeD &two) const { + return (x*two.x+y*two.y+z*two.z); +} + +template +void ThreeD::operator *= (const numtyp &two) { + x*=two; y*=two; z*=two; +} + +// Move coordinates into array +template +void ThreeD::to_array(numtyp *array) { + *array=x; + array++; + *array=y; + array++; + *array=z; +} + +// Set coordinates from array +template +void ThreeD::from_array(numtyp *array) { + x=*array; + array++; + y=*array; + array++; + z=*array; +} + + +// Return the cross product of *this and two +template +ThreeD ThreeD::cross(const ThreeD &two) const { + return ThreeD (y*two.z-z*two.y, z*two.x-x*two.z, x*two.y-y*two.x); +} + +// Return the angle between two vectors +template +numtyp ThreeD::angle(const ThreeD &two) const { + #ifdef NANCHECK + //assert(fabs(double((*this*two)/(norm()*two.norm())))<=1); + #endif + return(numtyp(acos( double(dot(two)/(norm()*two.norm())) ))); +} + +// Returns an arbitrary vector that is perpendicular to the input +template +ThreeD ThreeD::perpendicular() { + ThreeD two=*this; + if (y==0 && z==0) + two.y+=1; + else + two.x+=1; + return cross(two); +} + +template +ThreeD ThreeD::rotatey(double t) const { + ThreeD two; + numtyp sint=numtyp(sin(t)); + numtyp cost=numtyp(cos(t)); + two.x=(x*cost)-(z*sint); + two.z=(x*sint)+(z*cost); + two.y=y; + return two; +} + +template +ThreeD ThreeD::rotatez(double t) const { + ThreeD two; + numtyp sint=numtyp(sin(t)); + numtyp cost=numtyp(cos(t)); + two.x=(x*cost)-(y*sint); + two.y=(x*sint)+(y*cost); + two.z=z; + return two; +} + +template +ThreeD ThreeD::operator / (const numtyp &two) const { + return (ThreeD(x/two,y/two,z/two)); +} + +template +void ThreeD::operator /= (const numtyp &two) { + x/=two; y/=two; z/=two; +} + +// Magnitude of vector +template +numtyp ThreeD::hypot() const { + return norm(); +} + +// Return the norm of a vector +template +numtyp ThreeD::norm() const { + numtyp xi=x*x; numtyp yi=y*y; numtyp zi=z*z; + return numtyp(sqrt(double(xi+yi+zi))); +} + +// Return the squared norm of a vector +template +numtyp ThreeD::norm2() const { + return x*x+y*y+z*z; +} + +template +numtyp ThreeD::dist(const ThreeD &two) { + return (*this-two).norm(); +} + +template +numtyp ThreeD::dist2(const ThreeD &two) { + ThreeD diff=*this-two; + return diff.dot(diff); +} + +// For normalizing a vector +template +void ThreeD::normalize() { + numtyp temp=norm(); + #ifdef NANCHECK + assert(temp!=0); + #endif + *this/=temp; +} + +// Return unit vector +template +ThreeD ThreeD::unit() { + ThreeD unit=*this; + unit.normalize(); + return unit; +} + +/// Returns 3 +template +unsigned ThreeD::dimensionality() { + return 3; +} + +template +c2DPt ThreeD::projz() { + return c2DPt(x,y); +} + +template +ostream & operator<< (ostream &out, const ThreeD &t) { + out << t.x << " " << t.y << " " << t.z; + return out; +} + +template +istream & operator>>(istream &in, ThreeD &t) { + in >> t.x >> t.y >> t.z; + return in; +} + +template +ThreeD operator+ (const numtyp one, const ThreeD &two) { + return ThreeD(one+two.x,one+two.y,one+two.z); +} + +template +ThreeD operator- (const numtyp one, const ThreeD &two) { + return ThreeD(one-two.x,one-two.y,one-two.z); +} + +template +ThreeD operator* (const numtyp one, const ThreeD &two) { + return ThreeD(one*two.x,one*two.y,one*two.z); +} + +template +ThreeD operator/ (const numtyp one, const ThreeD &two) { + return ThreeD(one/two.x,one/two.y,one/two.z); +} + +// Set this to be the max of one and two for each dimension +template +void ThreeD::max(ThreeD &one, ThreeD &two) { + x=am::max(one.x,two.x); + y=am::max(one.y,two.y); + z=am::max(one.z,two.z); +} + +// Set this to be the max of one and two for each dimension +template +void ThreeD::min(ThreeD &one, ThreeD &two) { + x=am::min(one.x,two.x); + y=am::min(one.y,two.y); + z=am::min(one.z,two.z); +} + +// Returns the index of *this (for unsigned) in a square 3D array +/* s_size[0]=1Dsize, and s_size[1]=1D size*1Dsize */ +template +numtyp ThreeD::array_index(vector &s_size) { + return z+y*s_size[0]+x*s_size[1]; +} + +// Increment a 3D index from min to max (same as nested for) +template <> +bool ThreeD::increment_index(ThreeD &minp,ThreeD &maxp) { + x++; + if (x!=maxp.x) + return true; + x=minp.x; + y++; + if (y!=maxp.y) + return true; + y=minp.y; + z++; + if (z!=maxp.z) + return true; + return false; +} + +// Return false if x,y, or z is not within the inclusive range +template +bool ThreeD::check_bounds(numtyp minp,numtyp maxp) { + if (xmaxp || y>maxp || z>maxp) + return false; + return true; +} + +// --------------------- Quaternion Stuff ----------------------- + +Quaternion::Quaternion() { + #ifdef NANCHECK + w=0; i=0; j=0; k=0; + #endif +} + +Quaternion::~Quaternion() { +} + +Quaternion::Quaternion(const Quaternion &two) { + #ifdef NANCHECK + //assert(a::not_nan(two.w) && a::not_nan(two.i) && a::not_nan(two.j) && + // a::not_nan(two.k)); + #endif + w=two.w; i=two.i; j=two.j; k=two.k; +} + +Quaternion::Quaternion(double inw, double ini, double inj, double ink) { + #ifdef NANCHECK + //assert(a::not_nan(inw) && a::not_nan(ini) && a::not_nan(inj) && + // a::not_nan(ink)); + #endif + w=inw; i=ini; j=inj; k=ink; +} + +// Set the quaterion according to axis-angle (must be a UNIT vector) +Quaternion::Quaternion(const vectorPt &v, double angle) { + double halfa=angle/2; + double sina=sin(halfa); + w=cos(halfa); + i=v.x*sina; + j=v.y*sina; + k=v.z*sina; + #ifdef NANCHECK + //assert(a::not_nan(w) && a::not_nan(i) && a::not_nan(j) && a::not_nan(k)); + #endif +} + +// Set the quaterion according to rotation from vector 1 to vector 2 +// The input vectors do not need to be normalized +Quaternion::Quaternion(const vectorPt &v1, const vectorPt &v2) { + double angle=acos( v1.dot(v2)/(v1.norm()*v2.norm()) ); + vectorPt axis=v1.cross(v2); + double norm=axis.norm(); + if (norm>(istream &in, Quaternion &q) { + in >> q.w >> q.i >> q.j >> q.k; + return in; +} + +// ------------------- RotMat Stuff +RotMat::RotMat() { +} + +RotMat::RotMat(const Quaternion &q) { + set(q); +} + +void RotMat::set(const Quaternion &q) { + // x'=x(w^2+i^2-k^2-j^2)+y(2ij-2kw)+z(2jw+2ik) + // y'=x(2ij+2kw)+y(j^2-k^2+w^2-i^2)+z(2jk-2iw) + // z'=x(2ik-2jw)+y(2jk+2iw)+z(k^2-j^2-i^2+w^2) + double w2=q.w*q.w; + double i2=q.i*q.i; + double j2=q.j*q.j; + double k2=q.k*q.k; + double twoij=2.0*q.i*q.j; + double twoik=2.0*q.i*q.k; + double twojk=2.0*q.j*q.k; + double twoiw=2.0*q.i*q.w; + double twojw=2.0*q.j*q.w; + double twokw=2.0*q.k*q.w; + + x_x=w2+i2-j2-k2; + x_y=twoij-twokw; + x_z=twojw+twoik; + + y_x=twoij+twokw; + y_y=w2-i2+j2-k2; + y_z=twojk-twoiw; + + z_x=twoik-twojw; + z_y=twojk+twoiw; + z_z=w2-i2-j2+k2; +} + +cPt RotMat::operator*(const cPt &in) const { + cPt out; + out.x=x_x*in.x+x_y*in.y+x_z*in.z; + out.y=y_x*in.x+y_y*in.y+y_z*in.z; + out.z=z_x*in.x+z_y*in.y+z_z*in.z; + return out; +} + +void RotMat::proper() { + double det=x_x*y_y*z_z-x_x*y_z*z_y-x_y*y_x*z_z+x_y*y_z*z_x+x_z*y_x*z_y-x_z* + x_y*z_x; + if (det<0) { + x_x*=-1.0; + x_y*=-1.0; + x_z*=-1.0; + y_x*=-1.0; + y_y*=-1.0; + y_z*=-1.0; + z_x*=-1.0; + z_y*=-1.0; + z_z*=-1.0; + } +} + +EulerRot::EulerRot() { +} + +EulerRot::EulerRot(double theta_i,double psi_i,double phi_i) { + theta=theta_i; + psi=psi_i; + phi=phi_i; +} + +// Rotate the rotation axis by a given rotation matrix +void EulerRot::rotate_axis(const RotMat &rotmat) { + double c1=cos(theta/2.0); + double s1=sin(theta/2.0); + double c2=cos(psi/2.0); + double s2=sin(psi/2.0); + double c3=cos(phi/2.0); + double s3=sin(phi/2.0); + double c1c2=c1*c2; + double s1s2=s1*s2; + double w=c1c2*c3 - s1s2*s3; + + vectorPt axis; + axis.x=c1c2*s3 + s1s2*c3; + axis.y=s1*c2*c3 + c1*s2*s3; + axis.z=c1*s2*c3 - s1*c2*s3; + + double angle=2.0*acos(w); + double norm = axis.norm(); + if (norm < 0.00000001) + return; // No rotation + axis/=norm; + + axis=rotmat*axis; + + double s=sin(angle); + double c=cos(angle); + double t=1.0-c; + if ((axis.x*axis.y*t + axis.z*s) > 0.998) { // north pole singularity detected + theta = 2.0*atan2(axis.x*sin(angle/2.0),cos(angle/2.0)); + psi = PI/2.0; + phi = 0; + return; + } + if ((axis.x*axis.y*t + axis.z*s) < -0.998) { // south pole singularity detected + theta = -2.0*atan2(axis.x*sin(angle/2.0),cos(angle/2.0)); + psi = -PI/2.0; + phi = 0; + return; + } + + theta = atan2(axis.y*s-axis.x*axis.z*t, 1-(axis.y*axis.y+axis.z*axis.z)*t); + psi = asin(axis.x*axis.y*t+axis.z*s) ; + phi = atan2(axis.x*s-axis.y*axis.z*t, 1-(axis.x*axis.x+axis.z*axis.z)*t); +} + + +void EulerRot::operator= (const EulerRot &two) { + theta=two.theta; + psi=two.psi; + phi=two.phi; +} + +ostream & operator<<(ostream &out, const EulerRot &rot) { + out << rot.theta << " " << rot.psi << " " << rot.phi; + return out; +} + +// ------------------- Global functions + +// Find the point on the line in the direction of v closest to an outside point +// The point v_point is known to lie on the line +cPt c::point_to_line(const vectorPt &v, const cPt &v_point, const cPt &point){ + vectorPt u=point-v_point; + double uTv=u.dot(v); + double vTv=v.dot(v); + double qT=uTv/vTv; + + return ( (v*qT)+v_point); +} + +// Return the closest distance between two line segments input as end points +double c::closest_approach(const cPt &l1_1, const cPt &l1_2, const cPt &l2_1, + const cPt &l2_2) { + vectorPt u = l1_2 - l1_1; + vectorPt v = l2_2 - l2_1; + vectorPt w = l1_1-l2_1; + double a = u.dot(u); + double b = u.dot(v); + double c = v.dot(v); + double d = u.dot(w); + double e = v.dot(w); + double D = a*c - b*b; + double sc, sN, sD = D; + double tc, tN, tD = D; + + // compute the closest points between the two lines + if (D < 0.00000000001) { // parrallel lines + sN = 0.0; sD = 1.0; tN = e; tD = c; + } else { // get the closest points on the infinite lines + sN = (b*e - c*d); + tN = (a*e - b*d); + if (sN < 0.0) { + sN = 0.0; tN = e; tD = c; + } else if (sN > sD) { + sN = sD; tN = e + b; tD = c; + } + } + + if (tN < 0.0) { + tN = 0.0; + if (-d < 0.0) + sN = 0.0; + else if (-d > a) + sN = sD; + else { + sN = -d; sD = a; + } + } else if (tN > tD) { + tN = tD; + if ((-d + b) < 0.0) + sN = 0; + else if ((-d + b) > a) + sN = sD; + else { + sN = (-d + b); sD = a; + } + } + + sc=sN/sD; + tc=tN/tD; + + // vector connecting points + vectorPt dP=w+(u*sc)-(v*tc); + + return dP.norm(); +} + +void c::closest_approach_points(const cPt &l1_1, const cPt &l1_2, + const cPt &l2_1, const cPt &l2_2, + cPt &close_l1, cPt &close_l2) { + vectorPt u = l1_2 - l1_1; + vectorPt v = l2_2 - l2_1; + vectorPt w = l1_1-l2_1; + double a = u.dot(u); + double b = u.dot(v); + double c = v.dot(v); + double d = u.dot(w); + double e = v.dot(w); + double D = a*c - b*b; + double sc, sN, sD = D; + double tc, tN, tD = D; + + // compute the closest points between the two lines + if (D < 0.00000000001) { // parrallel lines + sN = 0.0; sD = 1.0; tN = e; tD = c; + } else { // get the closest points on the infinite lines + sN = (b*e - c*d); + tN = (a*e - b*d); + if (sN < 0.0) { + sN = 0.0; tN = e; tD = c; + } else if (sN > sD) { + sN = sD; tN = e + b; tD = c; + } + } + + if (tN < 0.0) { + tN = 0.0; + if (-d < 0.0) + sN = 0.0; + else if (-d > a) + sN = sD; + else { + sN = -d; sD = a; + } + } else if (tN > tD) { + tN = tD; + if ((-d + b) < 0.0) + sN = 0; + else if ((-d + b) > a) + sN = sD; + else { + sN = (-d + b); sD = a; + } + } + + sc=sN/sD; + tc=tN/tD; + + close_l1=u*sc+l1_1; + close_l2=v*tc+l2_1; + + return; +} + +// Returns true if two line segments intersect +bool c::intersect(const c2DPt &line1_start, const c2DPt &line1_end, + const c2DPt &line2_start, const c2DPt &line2_end) { + // Check segment1, line2 intersection + if (!sline_intersect(line1_start,line1_end,(line2_end-line2_start).normal(), + line2_start)) + return false; + + // Check segment2, line1 intersection + if (!sline_intersect(line2_start,line2_end,(line1_end-line1_start).normal(), + line1_start)) + return false; + + return true; +} + +// Liang-Barsky Intersection between line segment and line +bool c::sline_intersect(const c2DPt &line1_start, const c2DPt &line1_end, + const c2DPt &line2normal, const c2DPt &line2_point) { + double denom=(line1_end-line1_start).dot(line2normal); + if (fabs(denom)1) + return false; + return true; +} + +/// Average position +cPt c::mean(vector &vec) { + cPt average(vec[0]); + for (unsigned i=1; i +#include +#include +#include +using namespace std; + +enum { X, ///<0 + Y, ///<1 + Z ///<2 +}; +enum { I, ///<0 + J, ///<1 + K, ///<2 + W ///<3 +}; +enum { RED, ///<0 + GREEN, ///<1 + BLUE ///<2 +}; + +// Other coordinates +template class Ball; + +// Template friend declarations +template class TwoD; +template +ostream & operator<< (ostream &out, const TwoD &t); +template +istream & operator>> (istream &in, TwoD &t); + +template class ThreeD; +template +ostream & operator<< (ostream &out, const ThreeD &t); +template +istream & operator>> (istream &in, ThreeD &t); +template +ThreeD operator+ (const numtyp, const ThreeD &two); +template +ThreeD operator- (const numtyp, const ThreeD &two); +template +ThreeD operator* (const numtyp, const ThreeD &two); +template +ThreeD operator/ (const numtyp, const ThreeD &two); + +/// Two dimensional vector +/** The elements can be accessed directly .x or .y + * or by using the operator [] ( [X], [Y] or [I], [J] ) + * + * The following operators are currently overloaded: + * +,-,* + * + * operators *,/ returns a vector with each element multiplied + * times the corresponding element in the other vector (Matlab .* ) + * + * the member function dot can be used to compute dot products + * + * Input and output are overloaded for element I/O of the form "x y" + * <<, >> + * + * \sa cartesian.h for typedefs and defines*/ +template +class TwoD { + public: + /// Empty construct. Not necessarily initialized to [0 0] + TwoD(); + /// Assign both x and y the value + TwoD(numtyp x); + /// Assignment Constructor + TwoD(numtyp cx, numtyp cy); + + /// Type conversion + TwoD(const TwoD &two); + /// Ball projection (onto xy-plane) + TwoD(const Ball &ball); + + numtyp x; ///< First element + numtyp y; ///< Last element + + numtyp &operator[](unsigned i); + + friend ostream & operator<< <>(ostream &out, const TwoD &t); + friend istream & operator>> <>(istream &in, TwoD &t); + + TwoD operator + (const TwoD &two) const; + void operator += (const TwoD &two); + TwoD operator + (const numtyp two) const; + TwoD operator - (const TwoD &two) const; + TwoD operator - (const numtyp two) const; + TwoD operator * (const numtyp two) const; + TwoD operator * (const TwoD &two) const; + void operator /= (const numtyp two); + + bool operator != (const TwoD &two) const; + + /// Dot Product + numtyp dot(const TwoD &two) const; + + /// Distance between two points + numtyp dist(const TwoD &two) const; + /// Distance squared between two points + numtyp dist2(const TwoD &two) const; + /// Returns one of two normals to a line represented by vector + TwoD normal(); + + /// Move coordinates into array + void to_array(numtyp *array); + /// Set coordinates from array + void from_array(numtyp *array); + + // -------------- Weird functions that help with coord templating + unsigned dimensionality(); + /// Returns the index of *this (for unsigned) in a square 3D array + /** \param s_size s_size[0]=1Dsize **/ + numtyp array_index(vector &s_size); + /// Increment a 2D index from min to max (same as nested for) + /** Returns false when increment is complete **/ + bool increment_index(TwoD &minp,TwoD &maxp); + /// Return false if x or y is not within the inclusive range + bool check_bounds(numtyp min,numtyp max); + private: +}; + +///\var typedef TwoD c2DPt +/// Two dimensional vector of doubles +typedef TwoD c2DPt; + + +/// Three dimensional vector +/** The elements can be accessed directly .x or .y or .z + * or by using the operator [] ( [X], [Y], [Z] or [I], [J], [K] or + * [RED], [GREEN], [BLUE] ) + * + * The following operators are currently overloaded: + * +,-,*,/,+=,-=,*=,/=,==,!= + * + * operators *,/,*=,/= returns a vector with each element multiplied + * times the corresponding element in the other vector (Matlab .* ) + * or with each element multiplied by a scalar + * + * the member function dot can be used to compute dot products + + * + * Input and output are overloaded for element I/O of the form "x y z" + * <<, >> + * + * \sa cartesian.h for typedefs and defines*/ +template +class ThreeD { + public: + /// Assignment Constructor + ThreeD(numtyp cx, numtyp cy, numtyp cz); + /// Assign all the value + ThreeD(numtyp in); + /// Type Conversion + ThreeD(const ThreeD&); + /// Type Conversion + ThreeD(const ThreeD&); + /// Type Conversion + ThreeD(const ThreeD&); + /// Type Conversion + ThreeD(const ThreeD&); + /// TwoD with (z-coordinate set to zero) + ThreeD(const TwoD&); + /// Spherical Conversion + ThreeD(const Ball&); + /// Spherical Conversion + ThreeD(const Ball&); + /// Empty construct (Not necessarily initialized to zero) + ThreeD(); + + numtyp x; ///< First Element + numtyp y; ///< Second Element + numtyp z; ///< Last Element + + friend ostream & operator<< <>(ostream &out, const ThreeD &t); + friend istream & operator>> <>(istream &in, ThreeD &t); + friend ThreeD operator+ <>(const numtyp, const ThreeD &two); + friend ThreeD operator- <>(const numtyp, const ThreeD &two); + friend ThreeD operator* <>(const numtyp, const ThreeD &two); + friend ThreeD operator/ <>(const numtyp, const ThreeD &two); + + numtyp &operator[](unsigned i); + numtyp operator[](unsigned i) const; + + bool operator == (const ThreeD &two) const; + bool operator != (const ThreeD &two) const; + ThreeD operator + (const ThreeD &two) const; + ThreeD operator + (const numtyp &two) const; + ThreeD operator - (const ThreeD &two) const; + ThreeD operator - (const numtyp &two) const; + ThreeD operator * (const numtyp &two) const; + ThreeD operator * (const ThreeD &two) const; + ThreeD operator / (const numtyp &two) const; + ThreeD operator / (const ThreeD &two) const; + void operator = (const ThreeD &two); + void operator += (const numtyp &two); + void operator += (const ThreeD &two); + void operator -= (const numtyp &two); + void operator -= (const ThreeD &two); + void operator *= (const numtyp &two); + void operator /= (const numtyp &two); + + /// Move coordinates into array + void to_array(numtyp *array); + /// Set coordinates from array + void from_array(numtyp *array); + + /// Returns the dot product of *this and two + numtyp dot(const ThreeD &two) const; + + /// Returns the cross product of \b *this and \b two + /** The input vectors do not need to be normalized, however, the output + * vector will not be */ + ThreeD cross(const ThreeD &two) const; + /// Returns the angle between \b *this and \b two + numtyp angle(const ThreeD &two) const; + /// Returns an arbitrary vector that is perpendicular to the input + /** The output vector is not normalized */ + ThreeD perpendicular(); + /// Rotate a vector in xz-plane by t radians + ThreeD rotatey(double t) const; + /// Rotate a vector in xy-plane by t radians + ThreeD rotatez(double t) const; + /// Magnitude of vector + numtyp norm() const; + /// Squared norm of a vector + numtyp norm2() const; + /// Magnitude of vector + numtyp hypot() const; + /// Distance between two points + numtyp dist(const ThreeD &two); + /// Distance squared between two points + numtyp dist2(const ThreeD &two); + /// Converts \b *this to the unit vector + void normalize(); + /// Return the unit vector of \b *this + ThreeD unit(); + + /// Returns the projection of the point or vector onto Z-plane + c2DPt projz(); + + /// Set this to be the max of one and two for each dimension + void max(ThreeD &one, ThreeD &two); + /// Set this to be the min of one and two for each dimension + void min(ThreeD &one, ThreeD &two); + + // -------------- Weird functions that help with coord templating + /// Returns 3 + unsigned dimensionality(); + /// Returns the index of *this (for unsigned) in a square 3D array + /** \param s_size s_size[0]=1Dsize, and s_size[1]=1D size*1Dsize **/ + numtyp array_index(vector &s_size); + /// Increment a 3D index from min to max (same as nested for) + /** \note This is currently only implemented for unsigned numbers + * Returns false when increment is complete **/ + bool increment_index(ThreeD &minp,ThreeD &maxp); + /// Return false if x,y, or z is not within the inclusive range + bool check_bounds(numtyp min,numtyp max); + private: +}; + +///\var typedef ThreeD iPt; +/// Three dimensional vector of integers +typedef ThreeD iPt; +///\var typedef ThreeD uPt; +/// Three dimensional vector of unsigned +typedef ThreeD uPt; +///\var typedef ThreeD cPt; +/// Three dimensional vector of doubles +typedef ThreeD cPt; +///\var typedef ThreeD vectorPt; +/// Three dimensional vector of doubles +typedef ThreeD vectorPt; +///\var typedef ThreeD colorPt; +/// Three dimensional vector of doubles +typedef ThreeD colorPt; + +///\def ORIGIN +/// Point at origin +#define ORIGIN cPt(0.0,0.0,0.0) +///\def XAXIS +/// Unit vector for x-axis +#define XAXIS vectorPt(1.0,0.0,0.0) +///\def YAXIS +/// Unit vector for y-axis +#define YAXIS vectorPt(0.0,1.0,0.0) +///\def ZAXIS +/// Unit vector for z-axis +#define ZAXIS vectorPt(0.0,0.0,1.0) + +class RotMat; + +/// Euler Rotation +class EulerRot { + public: + EulerRot(); + EulerRot(double theta,double psi,double phi); + + /// Rotate the rotation axis by a given rotation matrix + void rotate_axis(const RotMat &rotmat); + + void operator= (const EulerRot &two); + friend ostream & operator<<(ostream &out, const EulerRot &rot); + + double theta; + double psi; + double phi; +}; + +///\def EU_NOROTATE +/// EulerRot with no rotation +#define EU_NOROTATE EulerRot(0.0,0.0,0.0) + +//---------------------------Quaternion Stuff ------------------------- + +/// Class for handling Quaternion vectors +/** The elements can be accessed directly .w .i .j or .k + * or by using the operator[] ( [W], [X], [Y], [Z] or [W], [I], [J], [K] ) + * + * The following operators are currently overloaded: + * +=,*,*=,=,== + * + * Overloaded * concatenates two successive rotations \n + * \e It \e is \e not \e communative \n + * For \e join=q1*q2, \e join is equivalent to performing rotation \e q1 + * \b followed by \e q2. + * + * Input and output are overloaded for element I/O of the form "w i j k" + * <<, >> + * + * \sa cartesian.h for typedefs and defines*/ +class Quaternion { + public: + Quaternion(); + ~Quaternion(); + + double w; ///>(istream &in, Quaternion &q); +}; + +/// Rotation matrix (about origin) +/** Rotation of ThreeD can be applied via overloaded * operator **/ +class RotMat { + public: + /// Unspecified matrix! + RotMat(); + /// Initialize with quaternion + RotMat(const Quaternion &q); + /// Set based on quaternion + void set(const Quaternion &q); + /// Assert that the rotation is proper + void proper(); + + cPt operator*(const cPt &in) const; + + double x_x,x_y,x_z,y_x,y_y,y_z,z_x,z_y,z_z; +}; + +///\def NOROTATE +/// Quaternion with no rotation +#define NOROTATE Quaternion(1.0,0.0,0.0,0.0) + +/// Global geometry functions +namespace c { + /// Returns point on a line closest to an outside point + /** The line is described by a directional vector and a point on the line + *\param v Vector describing the line direction + *\param v_point Point on the line + *\param point The Outside point */ + cPt point_to_line(const vectorPt &v, const cPt &v_point, + const cPt &point); + + /// Return the closest distance between two line segments input as end points + /**\param l1_1 End point of segment 1 + *\param l1_2 End point of segment 1 + *\param l2_1 End point of segment 2 + *\param l2_2 End point of segment 2 */ + double closest_approach(const cPt &l1_1, const cPt &l1_2, + const cPt &l2_1, const cPt &l2_2); + /// Calculates the points where two line segments are closest + /**\param l1_1 End point of segment 1 + *\param l1_2 End point of segment 1 + *\param l2_1 End point of segment 2 + *\param l2_2 End point of segment 2 + *\param close_l1 Calculated closest point on line segment 1 + *\param close_l2 Calculated closest point on line segment 2 */ + void closest_approach_points(const cPt &l1_1, const cPt &l1_2, + const cPt &l2_1, const cPt &l2_2, + cPt &close_l1, cPt &close_l2); + + /// Returns true if two line segments intersect + bool intersect(const c2DPt &line1_start, const c2DPt &line1_end, + const c2DPt &line2_start, const c2DPt &line2_end); + + /// Liang-Barsky Intersection between line segment and line + bool sline_intersect(const c2DPt &line1_start, const c2DPt &line1_end, + const c2DPt &line2normal, const c2DPt &line2_point); + /// Average position + cPt mean(vector &vec); +} + +#endif diff --git a/tools/pymol_asphere/src/colors.cpp b/tools/pymol_asphere/src/colors.cpp new file mode 100755 index 0000000000..4897442030 --- /dev/null +++ b/tools/pymol_asphere/src/colors.cpp @@ -0,0 +1,111 @@ +/*************************************************************************** + colors.cpp + W. Michael Brown + ------------------- + + begin : Wed Jun 11 2003 + copyright : (C) 2003 by mbrown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#include "colors.h" + +Colors::Colors() { + colors.insert(pair("white",colorPt(1,1,1))); + colors.insert(pair("red",colorPt(1,0,0))); + colors.insert(pair("green",colorPt(0,1,0))); + colors.insert(pair("blue",colorPt(0,0,1))); + colors.insert(pair("cmyk_blue",colorPt(0.074,0.168,0.614))); + colors.insert(pair("yellow",colorPt(1,1,0))); + colors.insert(pair("purple",colorPt(0.697,0.2,0.697))); + colors.insert(pair("hotpink",colorPt(1,0,0.5))); + colors.insert(pair("brown",colorPt(0.551,0.272,0.15))); + colors.insert(pair("orange",colorPt(1,0.5,0))); + colors.insert(pair("black",colorPt(0,0,0))); + colors.insert(pair("magenta",colorPt(1,0,1))); + colors.insert(pair("slate",colorPt(0.5,0.5,1))); + colors.insert(pair("marine",colorPt(0,0.5,1))); + colors.insert(pair("cmyk_marine",colorPt(0.273,0.469,0.758))); + colors.insert(pair("teal",colorPt(0.2,0,0.598))); + colors.insert(pair("forest",colorPt(0.098,0.5,0.098))); + colors.insert(pair("deep",colorPt(0.098,0.098,0.5))); + colors.insert(pair("grey",colorPt(0.5,0.5,0.5))); + colors.insert(pair("wheat",colorPt(0.988,0.819,0.646))); +} + +Colors::~Colors(){ +} + +// Number of colors in the map +unsigned Colors::size() { + return colors.size(); +} + +// Return the color in the gradient between start and end according to +// value +colorPt Colors::gradient(double start, double end, double value, + const colorPt &startcolor, const colorPt &endcolor) { + double percent; + + if (valueend) + return endcolor; + if (start==end) + return startcolor; + + percent=((value-start)/(end-start)); + return ((endcolor-startcolor)*percent)+startcolor; +} + +// Three color gradient +colorPt Colors::gradient(double start, double mid, double end, double value, + const colorPt &startcolor, const colorPt &midcolor, + const colorPt &endcolor) { + double percent; + + if (value==mid) + return midcolor; + else if (valueend) + return endcolor; + else if (value::iterator m; + + for (m=colors.begin(); m!=colors.end(); m++) + cout << " " << m->first << endl; +} + +// Return a string with the list of available colors +string Colors::colorlist() { + string colorl; + map::iterator m; + + for (m=colors.begin(); m!=colors.end(); m++) + colorl+="\n\t"+m->first; + colorl+="\n"; + return colorl; +} + +colorPt Colors::operator[](const string &name) { + map::iterator m; + m=colors.find(name); + if (m==colors.end()) + return colorPt(1,1,1); // Return white if not found + else + return m->second; +} + diff --git a/tools/pymol_asphere/src/colors.h b/tools/pymol_asphere/src/colors.h new file mode 100755 index 0000000000..af2fad1b6e --- /dev/null +++ b/tools/pymol_asphere/src/colors.h @@ -0,0 +1,59 @@ +/*************************************************************************** + colors.h + W. Michael Brown + ------------------- + + Storage and manipulation of colors. + + begin : Wed Jun 11 2003 + copyright : (C) 2003 by mbrown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#ifndef COLORS_H +#define COLORS_H + +#include "cartesian.h" +#include +#include +#include +using namespace std; + +/// Class for dealing with RGB color space +class Colors { +public: + Colors(); + ~Colors(); + + /// Return number of colors in the map + unsigned size(); + + /// Return the color in the gradient between start and end + /** \param start Starting value + * \param end Ending Value + * \param value Color in gradient determined for value + * relative to start and end + * \param startcolor Starting color for the gradient + * \param endcolor Ending color for the gradient */ + colorPt gradient(double start, double end, double value, + const colorPt &startcolor, const colorPt &endcolor); + /// Three color gradient + /** \sa gradient() */ + colorPt gradient(double start, double mid, double end, double value, + const colorPt &startcolor, const colorPt &midcolor, + const colorPt &endcolor); + + /// Output the colors to standard out + void outputcolors(); + /// Return a string with the list of available colors + string colorlist(); + + /// Return an RGB color based on name + colorPt operator[](const string &name); + +private: + map colors; +}; + +#endif + diff --git a/tools/pymol_asphere/src/commandline.cpp b/tools/pymol_asphere/src/commandline.cpp new file mode 100755 index 0000000000..cbf7bbc92f --- /dev/null +++ b/tools/pymol_asphere/src/commandline.cpp @@ -0,0 +1,504 @@ +/*************************************************************************** + commandline.cpp + W. Michael Brown + ------------------- + + Command line parsing stuff.. + * Add argument with ' ' for arguments that do not use - (i.e. filenames) + * Manditory for the ' ' requires only 1 argument regardless of the number + - use argsize; + + begin : Sun Jun 11 2003 + copyright : (C) 2003 by W. Michael Brown + email : mbrown@nirvana.unm.edu +***************************************************************************/ + +#include "commandline.h" + +CommandLine::Parameter::Parameter() { + num_args=0; + manditory_args=0; + manditory=false; + dash=false; + set=false; +} + +CommandLine::Parameter::~Parameter() {} + +CommandLine::CommandLine() { + help_set=false; + progname=""; + + // Set up the default man chapters + manchapter_order[0]="NAME"; + man_chapters["NAME"].clear(); + manchapter_order[1]="VERSION"; + man_chapters["VERSION"].clear(); + manchapter_order[2]="SYNOPSIS"; + man_chapters["SYNOPSIS"].clear(); + manchapter_order[3]="DESCRIPTION"; + man_chapters["DESCRIPTION"].clear(); + manchapter_order[4]="PARAMETERS"; + man_chapters["PARAMETERS"].clear(); + man_numchapters=5; +} + +CommandLine::~CommandLine() { +} + +void CommandLine::addargname(char n, const string &an) { + parameters[n].argnames.push_back(an); +} + +void CommandLine::adddescription(char n, const string &d) { + parameters[n].description.push_back(d); +} + +// Returns the argument size for a parameter +unsigned CommandLine::argsize(char n) { + return parameters[n].args.size(); +} + +bool CommandLine::operator [](char n) { + return set(n); +} + +string CommandLine::program_name() { + return progname; +} + +string CommandLine::full_command_line() { + return full_line; +} + +// ----------------- Add allowed arguments + +void CommandLine::add(char n, unsigned num) { + check(n); + parameters[n].num_args=num; + parameters[n].manditory_args=num; +} + +// Where num represents maximum number of arguments and man_num represents +// manditory number of arguments for a parameter +void CommandLine::add(char n, unsigned num, unsigned man_num) { + check(n); + parameters[n].num_args=num; + parameters[n].manditory_args=man_num; +} + +void CommandLine::addmanditory(char n, unsigned num) { + check(n); + parameters[n].manditory_args=num; + parameters[n].num_args=num; + parameters[n].manditory=true; +} + +// Where num represents maximum number of arguments and man_num represents +// manditory number of arguments for a parameter +void CommandLine::addmanditory(char n, unsigned num, unsigned man_num) { + check(n); + parameters[n].num_args=num; + parameters[n].manditory_args=man_num; + parameters[n].manditory=true; +} + +// Allow parameter that takes signed numbers as input +void CommandLine::addsigned(char n, unsigned num) { + check(n); + parameters[n].num_args=num; + parameters[n].manditory_args=num; + parameters[n].dash=true; // Allow dashes in arguments +} + +// Allow parameter that takes signed numbers as input +void CommandLine::addsigned(char n, unsigned num, unsigned man_num) { + check(n); + parameters[n].num_args=num; + parameters[n].manditory_args=man_num; + parameters[n].dash=true; // Allow dashes in arguments +} + +void CommandLine::addhelp(char n, unsigned num) { + check(n); + parameters[n].num_args=num; + parameters[n].manditory_args=num; + + help_set=true; + help_param=n; +} + +// Add descriptions for man pages +void CommandLine::addargnames(char n, unsigned num, const string args[]) { + for (unsigned i=0; i::iterator m; + m=parameters.find(n); + if (m!=parameters.end()) { + cerr << "DEVELOPER ERROR: Parameter: " << n << " set twice!\n"; + cerr << "commandline.h: Exiting...\n\n"; + exit(1); + } +} + +bool CommandLine::optargparams() { + map::iterator m; + for (m=parameters.begin(); m!=parameters.end(); m++) + if (m->second.manditory_argssecond.num_args) + return true; + + return false; +} + +// Returns the number of arguments that are not parameters (' ' name) +unsigned CommandLine::argsize() { + return parameters[' '].args.size(); +} + +// Whether or not this parameter was set +bool CommandLine::set(char n) { + return parameters[n].set; +} + +// Force a parameter to be unset +void CommandLine::unset(char n) { + parameters[n].set=false; +} + +char* CommandLine::arg(char n, unsigned num) { + return parameters[n].args[num]; +} + +int CommandLine::argint(char n, unsigned num) { + return atoi(parameters[n].args[num]); +} + +double CommandLine::argdouble(char n, unsigned num) { + return atof(parameters[n].args[num]); +} + +string CommandLine::argstring(char n, unsigned num) { + string s; + s=parameters[n].args[num]; + return s; +} + +bool CommandLine::parse(int argc, char* argv[], Error *error) { + unsigned i; + map::iterator m; + char flag; + int parameter; // Set to the parameter that arguments are being parsed + + progname=a::filenameonly(argv[0]); + full_line=string(argv[0]); + for (unsigned i=1; iaddwarning(3,9,"CommandLine","Invalid Command Line Argument: "+ + string(argv[argnum])); + return false; + } + parameters[' '].set=true; + parameters[' '].args.push_back(argv[argnum]); + argnum++; + continue; + } + + // Set a parameter + flag=argv[argnum][1]; + m=parameters.find(flag); + if (m==parameters.end()) { + error->addwarning(4,9,"CommandLine","Invalid Command Line Parameter: "+ + string(argv[argnum])); + return false; + } + // Make sure all required arguments are set before parameters with + // optional arguments + if (m->second.manditory_argssecond.num_args) + if (parameters[' '].args.size()buffer() << "Parameters with optional arguments must be " + << "placed after manditory commandline arguments"; + error->addbuf(5,9,"CommandLine"); + return false; + } + + parameter=argnum; + argnum++; + m->second.set=true; + + // Get the parameter arguments + for (i=0; isecond.num_args; i++) { + // Make sure we are not at the end of the parameter list + if (argnum>=argc) + if (m->second.args.size()second.manditory_args) { + error->buffer() << "Invalid Number of Arguments: -" << m->first + << " requires " << m->second.num_args + << " arguments!"; + error->addbuf(5,9,"CommandLine"); + return false; + } else + break; + + // Assure the right number of arguments for signed numbers + if (argv[argnum][0]=='-' && m->second.dash==true) + if (!isdigit(argv[argnum][1])) + if (m->second.args.size()second.manditory_args) { + error->buffer() << "Invalid Number of Arguments: -" << m->first + << " requires " << m->second.num_args + << " arguments!"; + error->addbuf(5,9,"CommandLine"); + return false; + } else + break; + + // Assure the right number of arguments for other numbers + if (argv[argnum][0]=='-' && m->second.dash==false) + if (m->second.args.size()second.manditory_args) { + error->buffer() << "Invalid Number of Arguments: -" << m->first + << " requires " << m->second.num_args + << " arguments!"; + error->addbuf(5,9,"CommandLine"); + return false; + } else + break; + + m->second.args.push_back(argv[argnum]); + argnum++; + } + } + + + // If help was set, we do not need to check for manditory args + if (help_set) + if (parameters[help_param].set) + return true; + + // Assure that the manditory arguments were set for commandline + if (parameters[' '].manditory) + if (parameters[' '].args.size()addwarning(6,9,"CommandLine","Missing Required Argument!"); + return false; + } + + // Assure that manditory arguments were set for parameters! + for (m=parameters.begin(); m!=parameters.end(); m++) { + if (m->second.manditory) + if (m->second.set==false) { + error->buffer() << "Missing Required Argument: \n\n"; + if (m->first!=' ') + error->buffer() << "-" << m->first << " must be set!"; + error->addbuf(7,9,"CommandLine"); + return false; + } + } + return true; +} + +void CommandLine::write_man_page(ostream & out, const string &version, + const string &header) { + unsigned i; + map::iterator m; + string bold="\\fB"; + string italic="\\fI"; + string regular="\\fR"; + + out << ".TH " << program_name() << " 1 \"" << a::date() + << "\"" << " \"" << program_name() << " (" << header << ") " + << version << "\" \"" << header << '\"' << endl; + + // Go through the chapters + map::iterator mc; + for (mc=manchapter_order.begin(); mc!=manchapter_order.end(); mc++) { + // NAME Section + if (mc->second=="NAME") { + out << ".SH NAME\n" << bold << program_name() << regular; + if (man_chapters["NAME"].size()!=0) { + out << " - "; + for (i=0; isecond=="SYNOPSIS" && man_chapters["SYNOPSIS"].empty()) { + out << ".PD 2\n.SH SYNOPSIS\n.PD 1\n.TP\n"; + if (program_name().length()>6) + out << bold << program_name() << regular << " "; + else + out << ".B " << program_name() << endl; + out << format_synopsis(bold,italic,regular) << endl; + if (program_name().length()>6) + out << ".br\n"; + if (parameters[' '].num_args!=0 && optargparams()) + out << ".PD 2\n.PP\nParameters with optional arguments should be placed " + << "after manditory commandline arguments.\n"; + continue; + } // end SYNOPSIS + + // PARAMETERS + if (mc->second=="PARAMETERS") { + if (!(parameters.size()==0 || + (parameters.size()==1 && (parameters.begin())->first==' '))) { + out << ".PD 2\n.SH PARAMETERS\n.PD 1\n"; + for (m=parameters.begin(); m!=parameters.end(); m++) { + if (m->first==' ') + continue; + out << ".TP\n"; + out << format_parameter(m->first,bold,italic,regular) << endl; + if (format_parameter(m->first,"","","").length()>6) + out << ".PD 0\n.TP\n.PP\n.PD 1\n"; + // Output the description + vector fdesc=man_format(m->second.description,bold, + italic,regular); + for (i=0; isecond,bold,italic,regular); + } +} + +void CommandLine::writeman_chapter(ostream &out,const string &name, + const string &bold, const string &italic, + const string ®ular) { + if (man_chapters[name].empty()) + return; + + out << ".PD 2\n.SH " << name << "\n.PD 1\n"; + + vector formatted=man_format(man_chapters[name],bold,italic,regular); + for (unsigned i=0; i::iterator m; + bool space=false; + + s=format_parameter(' ',bold,italic,regular)+" "; + for (m=parameters.begin(); m!=parameters.end(); m++) { + if (m->first==' ') + continue; + if (space) + s+=' '; + space=true; + if (m->second.manditory==false) + s+='['; + s+=format_parameter(m->first,bold,italic,regular); + if (m->second.manditory==false) + s+="]"; + } // end for + return s; +} + +// Write a synopsis in plain text format fitted to a given column width +void CommandLine::write_text_synopsis(ostream &out, unsigned column_width) { + string s=format_synopsis("","",""); + vector lines; + a::format_fit(column_width-program_name().length()-1,s,lines); + out << program_name() << " "; + for (unsigned i=0; i=p.manditory_args) + strp+='['; + if (param!=' ') + strp+=italic+p.argnames[i]+regular; + else + strp+=p.argnames[i]; + if (i>=p.manditory_args) + strp+=']'; + } + + if (p.argnames.size() CommandLine::man_format(const vector &input, + const string &bold, const string &italic, + const string ®ular) { + vector output; + map::iterator m; + + for (unsigned i=0; ifirst!=' ') { + string source="-"; + source+=m->first; + a::str_replace(source,bold+source+regular,temp); + } + for (unsigned j=0; jsecond.argnames.size(); j++) + a::str_replace(m->second.argnames[j], + italic+m->second.argnames[j]+regular,temp); + } + output.push_back(temp); + } + + return output; +} + diff --git a/tools/pymol_asphere/src/commandline.h b/tools/pymol_asphere/src/commandline.h new file mode 100755 index 0000000000..d0f4947218 --- /dev/null +++ b/tools/pymol_asphere/src/commandline.h @@ -0,0 +1,214 @@ +/*************************************************************************** + commandline.h + W. Michael Brown + ------------------- + + Command line parsing stuff.. + + begin : Sun Jun 11 2003 + copyright : (C) 2003 by W. Michael Brown + email : mbrown@nirvana.unm.edu + ***************************************************************************/ + +#ifndef COMMANDLINE_H +#define COMMANDLINE_H + +#include "error.h" +#include "misc.h" +#include +#include +#include +#include +#include +using namespace std; + +/// Parsing of command-line parameters and automatic man page generation +/** Allows manditory and optional command line arguments to be specified along + * with parameters set using flags. For arguments that do not require flags, + * a space char is used to set the flag. To specify the command line: + * + * foo input_file output_file -t signed_number [-o] + * + * \verbatim CommandLine cl; + cl.addmanditory(' ',2); + cl.addargname(' ',"input_file"); cl.addargname(' ',"output_file"); + cl.addsigned('t',1,1); + cl.addargname('t',"arg_name"); + cl.adddescription('t',"Manditory parameter allowing signed_numbers"); + cl.add('o',0); + \endverbatim + * To instead specify that the outputfile is an optional argument: + * \verbatim cl.addmanditory(' ',2,1) \endverbatim + * + * A help flag can also be set which does not require other command line + * arguments to be set + * + * When the commandline is parsed, it is asserted that all manditory arguments + * and flags are set, that each flag is set with the correct number of + * parameters, and that no unknown flags have been set. + * + * One can check whether or not optional flags have been set using set() or + * the [] operator: + * \verbatim cl.set('o') <--> cl['o'] \endverbatim + * + * Man pages can be generated automatically by adding chapters. The SYNOPSIS + * and description of commandline arguments and flags is generated + * automatically. Examples of chapters that are consistently used are + * NAME, VERSION, DESCRIPTION, PARAMETERS, USAGE, EXAMPLES, AUTHORS, BUGS, + * SEE ALSO. + * + * The following characters can be used for formatting in parameter + * descriptions and chapters: + * \verbatim + \t tab + \n forced newline + .TP new paragraph (margin indented) + \\fB all following text is bold + \\fI all following text is italic + \\fR all following text is regular + \endverbatim + * + * The arguments, flags and parameter_names are automatically formatted + * throughout the man page with bold and italic typeface. + **/ +class CommandLine { + public: + CommandLine(); + ~CommandLine(); + + /// Add an optional argument which takes num parameters + /** \note Use ' ' for n to specify an argument with no flag + * \param n flag used to pass arguments (i.e. -f) + * \param num the number of arguments that must be specified after flag */ + void add(char n, unsigned num); + /// Add a manditory argument which takes num parameters + /** \sa add() **/ + void addmanditory(char n, unsigned num); + /// Add a flag which can take signed numbers as parameters + /** \sa add() **/ + void addsigned(char n, unsigned num); + + /// Add an optional flag with minimum and maximum number of args + /** \sa add() **/ + void add(char n, unsigned num, unsigned man_num); + /// Add a manditory flag with minimum and maximum number of args + /** \sa add() **/ + void addmanditory(char n, unsigned num, unsigned man_num); + /// Add a flag with minimum and maximum number of args that takes signed nums + /** \sa add() **/ + void addsigned(char n, unsigned num, unsigned man_num); + + /// Add a help argument (does not require other manditory arguments be set) + /** \sa add() **/ + void addhelp(char n, unsigned num); + + /// Specify the names for arguments in the synopsis (for man page) + /** The names can be added one by one in the order of the flag parameters **/ + void addargname(char n, const string &an); + /// Specify the names for arguments in the synopsis (for man page) + void addargnames(char n, unsigned num, const string args[]); + /// Specify the description for an argument (for the man page SYNOPSIS) + /** See the class description for formating characters **/ + void adddescription(char n, const string &d); + /// Specify the description for an argument (for the man page SYNOPSIS) + /** See the class description for formating characters **/ + void adddescription(char n, unsigned num, const string d[]); + /// Add a man page chapter with title 'name' + void addtoman_chapter(const string &name,const string &body); + /// Add a man page chapter with title 'name' + void addtoman_chapter(const string &name,unsigned line_count, + const string body[]); + + /// Returns the number of arguments that are not flags (char n=' ') + unsigned argsize(); + /// Returns the number of parameters passed for a given flag + unsigned argsize(char n); + /// Returns true if the optional flag was set on the commandline + bool set(char n); + /// Returns true if the optional flag was set on the commandline + bool operator [](char n); + + /// Force a parameter to be unset + void unset(char n); + + /// Return flag parameter or argument as a cstring (0-based index) + char *arg(char n, unsigned num); + /// Return flag parameter or argument as a integer (0-based index) + int argint(char n, unsigned num); + /// Return flag parameter or argument as a double (0-based index) + double argdouble(char n, unsigned num); + /// Return flag parameter or argument as a string (0-based index) + string argstring(char n, unsigned num); + + /// Parse the command line arguments + bool parse(int argc, char* argv[], Error *error); + + /// Return the program name + string program_name(); + /// Return a string with the entire commandline as entered + string full_command_line(); + /// Write out a man_page + void write_man_page(ostream & out, const string &version, + const string &header); + + /// Advanced writing of man page + void writeman_chapter(ostream &out, const string &name, const string &bold, + const string &italic, const string ®ular); + + /// Return a string with all of the options in a man page synopsis format + string format_synopsis(const string &bold, const string &italic, + const string ®ular); + /// Write a synopsis in plain text format fitted to a given column width + void write_text_synopsis(ostream &out, unsigned column_width); + + private: + string full_line; + void check(char n); + // Returns true if parameters have been set with optional arguments + bool optargparams(); + + class Parameter; + map parameters; + + bool help_set; // True if there is a parameter for getting help + char help_param; // Parameter for help + + // Stuff for man page + string progname; // The name for the program (argv[0]) + unsigned man_numchapters; // Number of chapters in man page + map > man_chapters; + map manchapter_order; + + // Format a parameter as a string with argument names + string format_parameter(char param, const string &bold, const string &italic, + const string ®ular); + // Formats the strings in input for a man page by replacing newlines and + // adding bold and italic fonts to parameters and arguments respectively + vector man_format(const vector &input, const string &bold, + const string &italic, const string ®ular); + +}; + +// This class is for internal use within CommandLine +class CommandLine::Parameter { + public: + friend class CommandLine; + Parameter(); + ~Parameter(); + + private: + unsigned num_args; // Maximum number of arguments for the parameter + unsigned manditory_args; // Minimum number of arguments for the parameter + bool manditory; // This parameter MUST be set + + bool dash; // Allow for signed numbers (dashes in args) + + bool set; // This parameter has been set + vector args; + + // Strings for man page description + vector argnames; // Names for each argument to the parameter + vector description; // Description for how to set this parameter +}; + +#endif diff --git a/tools/pymol_asphere/src/error.cpp b/tools/pymol_asphere/src/error.cpp new file mode 100755 index 0000000000..ee6f4df0b7 --- /dev/null +++ b/tools/pymol_asphere/src/error.cpp @@ -0,0 +1,279 @@ +/*************************************************************************** + error.cpp + ------------------- + + Class for error handling + + __________________________________________________________________________ + + begin : Thu Oct 9 2003 + copyright : (C) 2003 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#include "error.h" + +Notice::Notice() { + nullout=new ostream(NULL); + noteout=&cout; + notice_level=9; +} + +Notice::~Notice() { + if (nullout!=NULL) + delete nullout; +} + +// Returns a null stream if level is two high, else returns notice stream +ostream & Notice::operator[] (const unsigned level) { + if (level<=notice_level) + return *noteout; + else + return *nullout; +} + +void Notice::setostream(ostream &out) { + noteout=&out; +} + +void Notice::set_notice_level(unsigned l) { + notice_level=l; +} + +unsigned Notice::get_notice_level() { + return notice_level; +} + +void Notice::notice(unsigned level, const string calling_class, + const string note) { + if (level ¬es) { + if (level<=notice_level) { + *noteout << calling_class << "::"; + for (unsigned i=0; i ¬es) { + if (level<=notice_level) + for (unsigned i=0; isecond.size(); +} + +void Error::addwarning(unsigned ID, unsigned level, const string calling_class, + const string warning) { + if (levelmax_level) + generate_error(ID,calling_class,warning); + vector *e=&(warning_list[ID]); + e->push_back(ErrCom()); + e->back().level=level; + e->back().calling_class=calling_class; + e->back().message=warning; + unhandled_warnings++; +} + +void Error::generate_error(unsigned ID, const string calling_class, + const string error) { + ErrCom err; + err.level=max_level+1; + err.calling_class=calling_class; + err.message=error; + + if (warnings()!=0) + writewarnings(); + write_err(ID,err); + writetotals(1); + #ifdef MUSE_MPI + MPI_Abort ( MPI_COMM_WORLD, 1 ); + #else + exit(1); + #endif +} + +// Add an error/warning (Program termination if level >= max level +ostringstream & Error::buffer() { + return buffer_stream; +} + +// Generate warning with message in buffer +void Error::addbuf(unsigned ID, unsigned level, const string calling_class) { + addwarning(ID,level,calling_class,buffer_stream.str()); + buffer_stream.str(""); +} + +// Generate serious error with message in buffer +void Error::addbuf(unsigned ID, const string calling_class) { + generate_error(ID,calling_class,buffer_stream.str()); +} + +unsigned Error::warnings() { + return unhandled_warnings; +} + +void Error::writeline() { + *errout << "+"; + *logout << "+"; + for (unsigned i=0; imax_level) { + et="| Error: "; width1=width12-7; + } else { + et="| Warning: "; width1=width12-9; + } + *errout << et << setw(width1) << ID << " | Level: " + << setw(width12-7) << err.level << " | " << setw(width3) + << err.calling_class << " |\n"; + *logout << et << setw(width1) << ID << " | Level: " + << setw(width12-7) << err.level << " | " << setw(width3) + << err.calling_class << " |\n"; + writeline(); + + // Output the message + vector messages; + a::format_fit(column_width-3,err.message,messages); + for (unsigned i=0; isecond.size()==1) + warning_list.erase(m); + else + m->second.erase(m->second.end()--); +} + +void Error::dismiss_all_warnings(unsigned ID) { + warning_iter m; + m=warning_list.find(ID); + if (m==warning_list.end()) + return; + unhandled_warnings-=m->second.size(); + warning_list.erase(m); +} + +void Error::writewarnings() { + while (warning_list.size()>0) + writewarning(warning_list.begin()->first); + return; +} + +void Error::dismiss_warnings() { + while (warning_list.size()>0) + dismiss_warning(warning_list.begin()->first); + return; +} + +// Write out the total warnings and errors +void Error::writetotals(unsigned errorcount) { + *errout << "\n"; + *logout << "\n"; + writeline(); + (*errout).setf(ios::left); + (*errout).unsetf(ios::right); + unsigned width1=(unsigned)floor((double)(column_width-7)/2.0); + unsigned width2=column_width-7-width1; + string swarnings="Total Warnings: "+a::itoa(handled_warnings+ + unhandled_warnings); + string serrors="Total Errors: "+a::itoa(errorcount); + *errout << "| " << setw(width1) << swarnings << " | " << setw(width2) + << serrors << " |\n"; + *logout << "| " << setw(width1) << swarnings << " | " << setw(width2) + << serrors << " |\n"; + writeline(); +} + diff --git a/tools/pymol_asphere/src/error.h b/tools/pymol_asphere/src/error.h new file mode 100755 index 0000000000..f14af6f253 --- /dev/null +++ b/tools/pymol_asphere/src/error.h @@ -0,0 +1,218 @@ +/*************************************************************************** + error.h + ------------------- + + Class for error handling + + __________________________________________________________________________ + + begin : Thu Oct 9 2003 + copyright : (C) 2003 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#ifndef ERRORCLASS +#define ERRORCLASS + +#include +#include +#include +#include +#include +#include +#include + +#ifdef MUSE_MPI +#include +#endif + +using namespace std; +// forward declarations +namespace a { + string itoa(unsigned int); + void format_fit(unsigned, const string &, vector &); +} + +/// Notice Class for Handling Object Output +/** A notice object stores an ostream for output and a notice_level. + * All output is sent along with a level. Any output whose level is + * greater than notice_level is sent to a null stream. The C++ output + * operator '<<' can be used with the Notice operator '[]' which passes + * the level: + * \verbatim notice_object[29] << "This notice has level 29" << endl; + * \endverbatim + * + * The guidelines for output notice levels are: + * - \e 0: Information that must be output + * - \e 1 - 9: Normal program output + * - \e 10-19: Parameters useful for storing how the program was run + * - \e 20-29: Extraneous information useful that may be useful to user + * - \e 30- : Debugging information */ +class Notice { + public: + /// Standard output (cout) is the default notice output ostream + /** The default maximum notice level is 9 \e (notice_level=10) */ + Notice(); + ~Notice(); + + /// Set the output stream for notice output + void setostream(ostream &out); + + /// Returns a null stream if level is two high, else returns notice stream + ostream & operator[] (const unsigned level); + + /// Set the degree of program output + void set_notice_level(unsigned l); + /// Get the degree of program output + unsigned get_notice_level(); + /// Generate a notice with a specified calling class + void notice(unsigned level, const string calling_class, const string note); + /// Generate a notice with a specified calling class + void notice(unsigned level, const string calling_class, + vector ¬es); + /// Generate a notice + void notice(unsigned level, const string note); + /// Generate a notice + void notice(unsigned level, vector ¬e); + + private: + unsigned notice_level; + + ostream *nullout; // Null stream for redirecting output to nowhere + ostream *noteout; // Output for notices +}; + +/// Error and Notice Handling +/** This class is intended to handle all output to the user. Output is + * divided into notices and warnings/errors. Output of any message is + * associated with a level. For notices, if level is greater than or equal to + * max_notice_level, no output is generated. For warnings, if level is less + * than min_warning_level, it is dismissed with no output. If the level is + * greater than max_warning_level, the program is terminated and all warnings + * and errors are output. + * + * \note By default, on destruction of an Error object, all unhandled + * warnings and errors are output + * + * A log file can be specified for each object. In this case, all notices + * are output to the log file only and errors are output to both stderr + * and the log file + * + * Errors can be generated with a string or using the internal message buffer: + \verbatim + Error error; + error.buffer() << "Incorrect file format for file: " << filename; + error.addbuf(512,19,"FooClass"; + // --- OR + string message = "Incorrect file format for file: "+filename; + error.addwarning(512,19,"FooClass",message); + \endverbatim + * + * Newlines will be inserted into the error message automatically in order + * to format it for the string. Forced newlines can be specified with \n + * + * Programs can check whether or not errors have been generated using the [] + * operator and can 'handle' them by outputing the message or dismissing + * them without any output + * + * Notices are generated using the public Notice class (see Notice()) + * + * \b Error \b Levels: + * - \e 0 - 1: Errors expected to happen during normal execution + * - \e 2 - 9: Errors that a non-interactive program can handle and continue + * - \e 10-19: Errors that interactive program can handle (file not found,etc.) + * - \e 20- : Serious errors that should terminate execution + **/ +class Error { + public: + /// Default constructor (use cerr for output and no log file) + /** Default max notice level is 9, min warning level is 2, and max warning + * level is 9 */ + Error(); + ~Error(); + + /// Set a log file for error AND notice output + void set_logfile(ostream &out); + + /// Returns the number of errors (if any) generated with id + unsigned operator[](unsigned id); + + /// Add warning, terminate if level is greater than max level + /** Newlines will be inserted into the message automatically when the + * message is formatted for output. However, forced newlines can also + * be inserted. **/ + void addwarning(unsigned ID, unsigned level, const string calling_class, + const string warning); + /// Add serious error (terminates execution) + void generate_error(unsigned ID, const string calling_class, + const string error); + + /// Add an message to the error buffer. Warning generated with addbuf() + /** Newlines will be inserted into the message automatically when the + * message is formatted for output. However, forced newlines can also + * be inserted. + * + \verbatim + Error error; + error.buffer() << "Choice not supported: " << choice; + error.addbuf(512,9,"FooClass"); + \endverbatim **/ + ostringstream & buffer(); + /// Generate warning with message in buffer + /** \sa buffer() **/ + void addbuf(unsigned ID, unsigned level, const string calling_class); + /// Generate serious error with message in buffer + /** \sa buffer() **/ + void addbuf(unsigned ID, const string calling_class); + + /// Number of Unhandled Warnings + unsigned warnings(); + /// Total number of warnings + unsigned total_warnings(); + + /// Handle all warnings with this ID by writing them to out + void writewarning(unsigned ID); + /// Handle LAST warning with this ID WITHOUT writing it + void dismiss_warning(unsigned ID); + /// Handle ALL warnings with this ID WITHOUT writing it + void dismiss_all_warnings(unsigned ID); + /// Handle all warnings by writing them out + void writewarnings(); + /// Handle all warnings without writing + void dismiss_warnings(); + + /// Write out the total warnings (write errorcount errors) + void writetotals(unsigned errorcount); + + /// For generating notices + /** \sa Notice **/ + Notice note; + private: + struct ErrCom; + map > warning_list; + typedef multimap >::iterator warning_iter; + unsigned handled_warnings; + unsigned unhandled_warnings; + bool handleatend; // Write any unhandled errors on destruct + bool writetotalatend; // Write totals on destruct if not 0 + + unsigned min_level; // Don't output warnings less than min_level + unsigned max_level; // if a warning has a level>max_level error! + ostream *errout, *logout; // Output for errors and warnings! + ostream *nullout; // No output + + ostringstream buffer_stream; // For creating messages for warnings + + unsigned column_width; + void write_err(unsigned ID, ErrCom &err); + void writeline(); +}; + +struct Error::ErrCom { + unsigned level; + string calling_class; + string message; +}; + +#endif + diff --git a/tools/pymol_asphere/src/glsurface.cpp b/tools/pymol_asphere/src/glsurface.cpp new file mode 100755 index 0000000000..23e484132a --- /dev/null +++ b/tools/pymol_asphere/src/glsurface.cpp @@ -0,0 +1,592 @@ +/*************************************************************************** + glsurface.cpp + W. Michael Brown + ------------------- + + Store and manipulate surfaces as OpenGL primitives + + begin : Sun Jun 8 2003 + copyright : (C) 2003 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#include "glsurface.h" + +GLSurface::GLSurface(){ +} +GLSurface::~GLSurface(){ +} + +// Operations on GLline structure +bool operator < (const GLline &one, const GLline &two) { + unsigned minone, mintwo, maxone, maxtwo; + + if (one.points[0] &t) { + triangles.swap(t); +} + +// Add many lines (via swap) +void GLSurface::addlinestrips(vector &l) { + linestrips.swap(l); +} + +// Add a mesh (via swap) +void GLSurface::addxyzmesh(vector &x,vector &y, + vector &z) { + xmesh.swap(x); ymesh.swap(y); zmesh.swap(z); +} + +// Preconditions: +// Triangles all loaded with correct vertex ordering +// Vertices all loaded +// Vertices that need to be calculated are marked as invalid +// Postconditions: +// For any normals that are marked as invalid, calculate normals by +// averaging normals for each triangle +void GLSurface::calculatenormals() { + unsigned i; + + // Set all invalid normals to zero + for (i=0; inote[31] << "Calculated normals: " << flush; + + return; +} + +// Normalize normals. +void GLSurface::normalize() { + unsigned i; + + for (i=0; inote[31] << "Normalized normals: " << flush; +} + +// Flip normals +void GLSurface::flipnormals() { + unsigned i; + + for (i=0; inote[31] << "Flipped normals: " << flush; +} + +// Calculate the Area of the Surface +double GLSurface::surfacearea() { + double sa=0; + + for (unsigned i=0; i &grid,const colorPt &neg, +// const colorPt &mid, const colorPt &pos, double minv, +// double midv, double maxv) { +// Colors c; +// for (unsigned i=0; i &l, ofstream &out, + const string &objname) { + unsigned i,j,k; + int index; + + writepymolheader(out); + //out << "LINEWIDTH, 2," << endl; + out << "BEGIN, LINES," << endl; + for (i=0; i +#include +#include +using namespace std; + +/// Stores the ID numbers of the vertices making a triangle +struct Triangle { + /// Triangle position + unsigned point[3]; +}; + +/// Vertex information used for all primitives +struct Vertex { + /// True if normal is being set for vertex + bool valid_normal; + /// Cartesian coords + cPt cpt; + // Used for surface area calculations + vectorPt normal; + /// Color + colorPt color; + /// For coloring surface + char marker; + /// For transparency (0 transparent, 1 opaque) + double transparency; +}; + +/// ID numbers for vertices making a line +struct GLline { + unsigned points[2]; +}; + +// This stuff is used for removing duplicate lines by sorting; +bool operator < (const GLline &one, const GLline &two); +bool operator == (const GLline &one, const GLline &two); + +/// Line Strip (loop) primitive (stores vertex indices) +struct LineStrip { + /// True for loop, false for line_strip + bool loop; + vector line; +}; + +/// Sphere primitive +struct Sphere { + unsigned i; + double radius; +}; + +/// Class for storage and I/O of a set of OpenGL primitives +class GLSurface { +public: + GLSurface(); + ~GLSurface(); + /// Reserve vector space for specified number of vertices and triangles + void reserve(unsigned num_v, unsigned num_t); + /// Delete all primitives + void clear(); + + /// Add a triangle + void addtriangle(Triangle &t); + /// Add a vertex + void addvertex(Vertex &v); + + /// Add many triangles (via swap) + /** \note triangles are removed from input vector */ + void addtriangles(vector &t); + /// Add many lines (via swap) + /** \note lines are removed from input vector */ + void addlinestrips(vector &l); + /// Add a mesh (via swap) + /** \note lines are removed from input vector */ + void addxyzmesh(vector &x,vector &y,vector &z); + /// Add a line + void addline(unsigned v1, unsigned v2); + /// Add a sphere based on a current vertex and a radius + void add_sphere(unsigned index, double radius) { + Sphere s; + s.i=index; + s.radius=radius; + spheres.push_back(s); + } + + /// Add an ellipsoid with specified resolution + /** The number of triangles is equal to twice the resolution squared **/ + void add_ellipsoid(const cPt ¢er, const cPt &rad, + const Quaternion &rot, const colorPt &color, + const double alpha,const unsigned resolution); + /// Add an ellipsoid with resolution of 10 + /** \sa add_ellipsoid **/ + void add_ellipsoid(const cPt ¢er, const cPt &rad, + const Quaternion &rot, const colorPt &color, + const double alpha); + + /// Calculate the unit normals for any vertices marked with invalid normals + void calculatenormals(); + /// Normalize normals. + void normalize(); + /// Flip normals + void flipnormals(); + /// Calculate the Area of the Surface + double surfacearea(); + /// Calculate the Area of a Triangle + double triarea(Triangle &t); + + /// Return the number of spheres + unsigned size_spheres() { return spheres.size(); } + /// Return the number of vertices + unsigned size_vertices() { return vertices.size(); } + + /// Solid color a surface + void color(const colorPt &color); + /// Color a surface by using a gradient + void colorbygradient(const colorPt &start, const colorPt &end); + /// Color a surface based on a marker + void colorbymarker(const colorPt &zero, const colorPt &nonzero); + /// Color a surface based on interpolation of values on a grid + /** \param grid Grid with appropriate interpolation scheme set + * \param neg Color for the value at and below minv + * \param mid Color for the value at midv + * \param pos Color for the value at and above posv **/ + //void colorbygrid(GridNUM &grid,const colorPt &neg, + // const colorPt &mid, const colorPt &pos, double minv, + // double midv, double maxv); + /// Set the transparency of the entire surface (0-1 [1=no transparency]) + void set_transparency(double alpha); + + /// Write out triangles as BEGIN,TRIANGLES,END primitives + void writetris(ofstream &out, const string &objname); + /// Write out triangles as TRIANGLE primitives + void writetris_surf(ofstream &out, const string &objname); + /// Write out triangles as POINTS + void writetris_points(ofstream &out, const string &objname); + /// Write out vertices as spheres + void write_vspheres(ofstream &out, const string &objname, double radius); + /// Write out triangles as Mesh + void writetris_mesh(ofstream &out, const string &objname); + + /// Write out triangles as an xyzmesh + void writexyzmesh(ofstream &out,const string &objname); + + void writelines(ofstream &out, const string &objname); + void writelines(vector &l, ofstream &out, const string &objname); + /// Write out Line Strips (Loops) + void writelinestrips(ofstream &out, const string &objname); + void writespheres(ofstream &out, const string &objname); + + /// Write header for Python scripts for rendering in PyMol + void writepymolheader(ofstream &out); + /// Write tail for Python script for rendering in Pymol + void writepymoltail(ofstream &out, const string &objname); +private: + vector triangles; + vector vertices; + vector gllines; + vector spheres; + + // For xyzmesh + vector xmesh; + vector ymesh; + vector zmesh; + + vector linestrips; + + void add_super_ellipsoid(const cPt &cen, const cPt &rad,const Quaternion &q, + const double n, const double e, const double u1, + const double u2, const double v1, const double v2, + const unsigned u_segs, const unsigned v_segs, + const colorPt &color, const double alpha); + void SQE_helper(Vertex &ver, const cPt &rad, const double u, + const double v, const double n, const double e); +}; + +#endif + diff --git a/tools/pymol_asphere/src/m_constants.h b/tools/pymol_asphere/src/m_constants.h new file mode 100755 index 0000000000..ca08206526 --- /dev/null +++ b/tools/pymol_asphere/src/m_constants.h @@ -0,0 +1,67 @@ +/*************************************************************************** + m_constants.h + ------------------- + W. Michael Brown + + Misc constants + + __________________________________________________________________________ + This file is part of the Math Library + __________________________________________________________________________ + + begin : Wed Aug 10 2005 + copyright : (C) 2005 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +/*! \file */ + +#ifndef CONSTANTS_H +#define CONSTANTS_H + +#include +using namespace std; + +#define MATHLIB_VER "0.15" + +#ifndef PI +///\def PI +/// pi +#define PI 3.1415926535897932384626433832795 +#endif +///\def TWOPI +/// pi*2 +#define TWOPI 6.2831853071795862 +///\def HALFPI +/// pi/2 +#define HALFPI 1.5707963267948966 +///\def DEGTORAD +/// Convert Degrees to Radians (pi/180) +#define DEGTORAD 0.017453292519943295 +///\def SQRT_TWO +/// sqrt(2.0) +#define SQRT_TWO 1.4142135623730951 +///\def SQRT_PI +/// sqrt(PI) +#define SQRT_PI 1.7724538509055159 +///\def INF +/// Infinity +#define INF 1e308 +///\def MINUSINF +/// Negative infinity +#define MINUSINF -1e308 + +#ifndef EPS +///\def EPS +/// Small number +#define EPS 1e-100 +#endif + +/** \mainpage Math Library + * \section intro Introduction + * Math library with containers and operations for vectors, matrices, graphs, + * cartesian coordinates, quaternions, Euler angles, support vector machine + * models, etc. **/ + + +#endif diff --git a/tools/pymol_asphere/src/make_manpages.sh b/tools/pymol_asphere/src/make_manpages.sh new file mode 100755 index 0000000000..678d059b4f --- /dev/null +++ b/tools/pymol_asphere/src/make_manpages.sh @@ -0,0 +1,24 @@ +#!/bin/tcsh + +set execs="asphere_vis" +set execdir="../bin" + +if ( -e ../doc ) then + echo Manpage directory exists... +else + echo Creating directory 'manpages' + mkdir ../doc +endif + +cd ../doc + +foreach exec ($execs) + $execdir/$exec -h > $exec.manpage + eqn $exec.manpage > $exec.1 + man -t -p eqn ./$exec.manpage > $exec.ps + ps2pdf $exec.ps $exec.pdf +end + +cd ../src + + diff --git a/tools/pymol_asphere/src/misc.cpp b/tools/pymol_asphere/src/misc.cpp new file mode 100755 index 0000000000..8696f4410d --- /dev/null +++ b/tools/pymol_asphere/src/misc.cpp @@ -0,0 +1,428 @@ +/*************************************************************************** + misc.cpp + ------------------- + W. Michael Brown + + Miscellaneous functions that do not deserve their own class + + __________________________________________________________________________ + This file is part of the "All" Library + __________________________________________________________________________ + + begin : May 30 2003 + copyright : (C) 2003 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#include "misc.h" + +void a::copyright(ostream &out, unsigned year) { + out << "Copyright (" << year << ") Sandia Corporation. Under the terms of " + << "Contract\nDE-AC04-94AL85000, there is a non-exclusive license for " + << "use of this\nwork by or on behalf of the U.S. Government. Export " + << "of this program\nmay require a license from the United States " + << "Government.\n"; +} + +// Get the SNL Copyright info as a string +string a::copyrightstring(unsigned year) { + return string("Copyright (")+a::itoa(year)+") Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive license for use of this work by or on behalf of the U.S. Government. Export of this program may require a license from the United States Government."; +} + +void a::fileopen(ifstream &in, const char *filename, Error &error) { + in.clear(); + in.open(filename); + if (!in) { + error.addwarning(1,15,"Misc", + "Could not open "+string(filename)+" for input!"); + error.writewarnings(); + } +} + +void a::fileopen(ifstream &in, const string &filename, Error &error) { + in.clear(); + in.open(filename.c_str()); + if (!in) { + error.addwarning(1,15,"Misc", + "Could not open "+string(filename)+" for input!"); + error.writewarnings(); + } +} + +void a::fileopenbinary(ifstream &in, const string &filename, Error &error) { + in.clear(); + in.open(filename.c_str(), ios::binary); + if (!in) { + error.addwarning(1,15,"Misc", + "Could not open "+string(filename)+" for input!"); + error.writewarnings(); + } +} + +void a::fileopen(ofstream &out, const char *filename, Error &error) { + out.clear(); + out.open(filename); + if (!out) { + error.addwarning(2,15,"Misc", + "Could not open "+string(filename)+" for output!"); + error.writewarnings(); + } +} + +void a::fileopen(ofstream &out, const string &filename, Error &error) { + out.clear(); + out.open(filename.c_str()); + if (!out) { + error.addwarning(2,15,"Misc", + "Could not open "+string(filename)+" for output!"); + error.writewarnings(); + } +} + +void a::fileopenbinary(ofstream &out, const string &filename, Error &error) { + out.clear(); + out.open(filename.c_str(),ios::binary); + if (!out) { + error.addwarning(2,15,"Misc", + "Could not open "+string(filename)+" for output!"); + error.writewarnings(); + } +} + +void a::fileopenapp(ofstream &out, const string &filename, Error &error) { + out.clear(); + out.open(filename.c_str(),ios::app); + if (!out) { + error.addwarning(2,15,"Misc", + "Could not open "+string(filename)+" for output!"); + error.writewarnings(); + } +} + +void a::fileclose(ifstream &in, Error &error) { + in.close(); +} + +void a::fileclose(ofstream &out, Error &error) { + if (out.fail()) { + error.addwarning(10,15,"Misc","Error writing to output file!"); + error.writewarnings(); + } + out.close(); +} + +// Put a string back into an istream +void a::putback(istream &in, const string &s) { + if (s.size()==0) + return; + unsigned i=s.size()-1; + while (true) { + in.putback(s[i]); + if (i==0) + return; + i--; + } +} + +string a::date() { + char datestr[40]; + + time_t t; + time(&t); + + strftime(datestr,40,"%B %d, %Y",localtime(&t)); + return string(datestr); +} + +// Return the filename without the extension +string a::namewoext(const string &filename) { + return (filename.substr(0,filename.find_last_of('.'))); +} + +// Return the filename without extension or directory +string a::filenameonly(const string &filename) { + // Find the start of the filename + unsigned start=0; + if (filename.find_last_of('/')<(filename.length()-1)) + start=filename.find_last_of('/')+1; + + // Find the end of the filename + unsigned end=filename.find_last_of('.')-start; + return(filename.substr(start,end)); +} + +// Return the extension of a filename +string a::extension(const string &filename) { + // Find the start of the extension + unsigned start=filename.find_last_of('.')+1; + if (start>=filename.length()) + return ""; + + return filename.substr(start,filename.length()-start); +} + +// Center a string over the specified length +string a::strcenter(const string &s, unsigned length) { + string empty(" "); + unsigned half=length/2; + unsigned spacer=half-(s.length()/2)-1; + return (empty.substr(0,spacer)+s); +} + +// True if a character is whitespace +bool a::whitespace(char c) { + if (c==' ' || c=='\n' || c=='\t' || c=='\f' || c=='\r' || c=='\v') + return true; + return false; +} + +// Check if a string is only whitespace +bool a::whitespace(const string &s) { + for (unsigned i=0; i0 && tlength>0); + #endif + + while (true) { + if (loc>=s.length()) + break; + loc=s.find(source,loc); + if (loc>=s.length()) + break; + s.replace(loc,slength,target,0,tlength); + loc+=tlength; + } +} + +/// Convert all alpha characters to lower case +string a::tolower(const string &s) { + string r=""; + for (unsigned i=0; i &tokens) { + string sline=line; + get_tokens(sline,tokens); +} + +void a::get_tokens(const string &sline, vector &tokens) { + string token=""; + unsigned i=0; + while (i &tokens) { + string token=""; + unsigned i=0; + while (i &output) { + vector forced; + a::get_tokens('\n',input,forced); + for (unsigned i=0; i tokens; + a::get_tokens(forced[i],tokens); + for (unsigned j=0; jcolumn_width) { + unsigned this_count=column_width-current_line.length(); + current_line+=tokens[j].substr(0,this_count); + output.push_back(current_line); + current_line=""; + tokens[j]=tokens[j].substr(this_count,tokens[j].length()-this_count); + j--; + } else { + output.push_back(current_line); + current_line=tokens[j]+' '; + } + } + } + output.push_back(current_line); + } + return; +} + +string a::itoa(unsigned i) { + ostringstream o; + o << i; + return o.str(); +} + +string a::itoa(int i) { + ostringstream o; + o << i; + return o.str(); +} + +string a::ftoa(double i) { + ostringstream o; + o << i; + return o.str(); +} + +// Seed the random number generator +void a::seedrandom(unsigned seed) { + srand(seed); +} + +// Seed the random number generator with the current time +void a::seedrandom_time() { + srand(unsigned(time(0))); //+getpid()); +} + +// Return an integer between 0 and max +double a::frandom(double max) { + return double(rand())/double(RAND_MAX)*max; +} + +// Return an integer between 0 and max +unsigned a::irandom(unsigned max) { + return unsigned(double(rand())/double(RAND_MAX)*max); +} + +// Return an integer between 0 and max +long a::lrandom(long max) { + return long(double(rand())/double(RAND_MAX)*max); +} + +// Empty constructer with no header, no extension, no lead zeros +FileIterator::FileIterator() { + file_num=0; + digits=0; + header=""; + extension=""; +} + +// Specify the filename format with leading zeros +FileIterator::FileIterator(const string &h,const string &e,unsigned d) { + file_num=0; + digits=d; + header=h; + extension=e; +} + +// Specify the filename format without leading zeros +FileIterator::FileIterator(const string &h,const string &e) { + digits=0; + file_num=0; + header=h; + extension=e; +} + +// Set the current file number +void FileIterator::set_file_num(unsigned fnum) { + file_num=fnum; +} + +// Set the file header +void FileIterator::set_file_header(const string &head) { + header=head; +} + +// Set the file extension +void FileIterator::set_file_extensions(const string &ext) { + extension=ext; +} + +// Set the number of leading zeros +void FileIterator::set_lead_zeros(unsigned digs) { + digits=digs; +} + +// Set the current file number, header, extension, leading zeros +void FileIterator::set(unsigned fnum, const string &head, const string &ext, + unsigned digs) { + file_num=fnum; + header=head; + extension=ext; + digits=digs; +} + +string FileIterator::nextfilename() { + // Get an output filename from a string + string filename; + filename=a::itoa(file_num); + while (filename.length() +#include +#include +#include +#include +#include +#include + +using namespace std; + +/// Miscellaneous functions that do not deserve their own class +/** \e a contains functions for \n + * - fileio + * - string manipulation + * - conversion from numbers to strings */ +namespace a { + /// Output the SNL Copyright info for publication \e year. + void copyright(ostream &out, unsigned year); + /// Get the SNL Copyright info as a string + string copyrightstring(unsigned year); + + /// Open a file for input. Generates error ID \b 1-15 on fail. + void fileopen(ifstream &in, const char *filename, Error &error); + /// Open a file for input. Generates error ID \b 1-15 on fail. + void fileopen(ifstream &in, const string &filename, Error &error); + /// Open a binary file for input. Generates error ID \b 1-15 on fail. + void fileopenbinary(ifstream &in, const string &filename, Error &error); + + /// Open a file for output. Generates error ID \b 2-15 on fail. + void fileopen(ofstream &out, const string &filename, Error &error); + /// Open a file for output. Generates error ID \b 2-15 on fail. + void fileopenbinary(ofstream &out, const string &filename, Error &error); + /// Open a file for append. Generates error ID \b 2-15 on fail. + void fileopenapp(ofstream &out, const string &filename, Error &error); + /// Open a file for output. Generates error ID \b 2-15 on fail. + void fileopen(ofstream &out, const char *filename, Error &error); + + /// Close an input file. Generates error ID \b 10-15 on fail. + void fileclose(ifstream &in, Error &error); + /// Close an output file. Generates error ID \b 11-15 on fail. + void fileclose(ofstream &out, Error &error); + + /// Put a string back into an istream + void putback(istream &in, const string &s); + + /// Get the current date in the format "January 1, 2003" + string date(); + + /// Returns the filename without the extension + string namewoext(const string &filename); + /// Returns the filename without extension or directory + string filenameonly(const string &filename); + /// Return the extension of a filename + string extension(const string &filename); + /// Centers a string over the specified length + string strcenter(const string &s, unsigned length); + /// True if a character is whitespace + bool whitespace(char c); + /// True if a string is only whitespace + bool whitespace(const string &s); + /// Remove all whitespace from a string + string remove_whitespace(const string &s); + /// Replace any instance of \e source with \e target within the string \e s + void str_replace(const string &source, const string &target, string &s); + /// Convert all alpha characters to lower case + string tolower(const string &s); + + /// The tokens parsed from cstring are \e added to the input vector + /** \param line The line with 'white space' delimeted tokens + * \param tokens Each token parsed is added to the vector */ + void get_tokens(const char *line, vector &tokens); + /// The tokens parsed from string are \e added to the input vector + /** \param line The line with 'white space' delimeted tokens + * \param tokens Each token parsed is added to the vector */ + void get_tokens(const string &line, vector &tokens); + /// Parse a string into tokens based on delimiter + /** \param line The line with 'white space' delimeted tokens + * \param tokens Each token parsed is added to the vector */ + void get_tokens(char delimiter, const string &line, vector &tokens); + /// Return the first token in a string + string get_first_token(const char *line); + /// Format a string to fit within a specified column width + /** Newlines are inserted between whitespace if possible, otherwise line is + * wrapped. Each string in the vector represents one line of text + * \note The output vector is not emptied! **/ + void format_fit(unsigned column_width, const string &input, + vector &output); + + /// Return a string of num underscores + string underline(unsigned num); + + /// Returns string representation of unsigned number + string itoa(unsigned i); + /// Returns string representation of int number + string itoa(int i); + /// Returns string representation of double number + string ftoa(double i); + + /// Seed the random number generator + void seedrandom(unsigned seed); + /// Seed the random number generator with the current time + void seedrandom_time(); + /// Returns a random integer between 0 and max + unsigned irandom(unsigned max); + /// Returns a random integer between 0 and max + long lrandom(long max); + /// Returns a random double between 0 and max + double frandom(double max); +} + +/// Iterate through file names to give each unique numbers +class FileIterator { + public: + /// Empty constructer with no header, no extension, no lead zeros + FileIterator(); + /// Specify the filename format using leading zeros + /** Files are generated according to the following format: + * \verbatim header+%0'digits'file_number+extension \endverbatim */ + FileIterator(const string &header, const string &extension, unsigned digits); + /// Specify the filename format without leading zeros + /** Files are generated according to the following format: + * \verbatim header+file_number+extension \endverbatim */ + FileIterator(const string &header, const string &extension); + + /// Set the current file number + void set_file_num(unsigned fnum); + /// Set the file header + void set_file_header(const string &head); + /// Set the file extension + void set_file_extensions(const string &ext); + /// Set the number of leading zeros + void set_lead_zeros(unsigned digs); + /// Set the current file number, header, extension, leading zeros + void set(unsigned fnum, const string &head, const string &ext, unsigned digs); + + /// Returns the next filename. + string nextfilename(); + private: + string header,extension; + unsigned digits; + unsigned file_num; +}; + +#endif diff --git a/tools/pymol_asphere/src/miscm.cpp b/tools/pymol_asphere/src/miscm.cpp new file mode 100755 index 0000000000..61174da67c --- /dev/null +++ b/tools/pymol_asphere/src/miscm.cpp @@ -0,0 +1,212 @@ +/*************************************************************************** + miscm.cpp + ------------------- + W. Michael Brown + + Miscellaneous functions that do not deserve their own class + + __________________________________________________________________________ + This file is part of the Math Library + __________________________________________________________________________ + + begin : May 30 2003 + copyright : (C) 2003 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#include "miscm.h" + +double am::square(double num) { + return num*num; +} + +// Rounds a number +double am::round(double n) { + double r=ceil(n); + if (r-n>0.5) + return floor(n); + return r; +} + +// Return the -1 for negative, 0 for zero, and 1 for positive +int am::sign(double v) { + if (v<0) + return -1; + if (v>0) + return 1; + return 0; +} + +// Return the range of elements in a vector +double am::range(const vector &v) { + if (v.empty()) + return 0; + double min=v[0]; + double max=v[0]; + for (unsigned i=1; imax) + max=v[i]; + } + return max-min; +} + +// Return the average of elements in a vector +double am::mean(const vector &v) { + double sum=0; + for (unsigned i=0; i(double,double); + template float max(float,float); + template unsigned max(unsigned,unsigned); + template int max(int,int); +} + +template +numtyp am::max(numtyp one,numtyp two) { + if (one>two) + return one; + return two; +} + +// Return the min of two objects +namespace am { + template double min(double,double); + template float min(float,float); + template unsigned min(unsigned,unsigned); + template int min(int,int); +} + +template +numtyp am::min(numtyp one,numtyp two) { + if (one= big) { + big =fabs(a[j][k]); + irow=j; + icol=k; + } + } else if (ipiv[k] > 1) { + error.addwarning(303,9,"Invert", + "Cannot invert a singular matrix."); + return; + } + } + ++(ipiv[icol]); + if (irow !=icol) { + for (l=0;l=1;l--) { + if (indxr[l] != indxc[l]) + for (k=0;k +#include +#include +#include "error.h" + +using namespace std; + +/// Miscellaneous functions that do not deserve their own class +/** \e a contains functions for \n + * - simple math functions + * - finite precision stuff */ + +namespace am { + /// Returns the square of a number + double square(double); + /// Rounds a number + double round(double); + + /// Return the range of elements in a vector (0 if empty) + double range(const vector &v); + /// Return the average of elements in a vector (0 if empty) + double mean(const vector &v); + + /// Return the max of two objects + template + numtyp max(numtyp one,numtyp two); + + /// Return the min of two objects + template + numtyp min(numtyp one,numtyp two); + + /// Return the -1 for negative, 0 for zero, and 1 for positive + int sign(double v); + + /// Swap two objects + void swap(double a, double b); + + // Finite Precision stuff + + /// Returns a number representing zero for finite checks + double epsilon(double number); + /// Returns number closer to zero by the smallest interval possible + double minus_eps(double number); + /// Returns number farther from zero by the smallest interval possible + double plus_eps(double number); + /// Returns number farther from zero by smallest interval * \b m + double plus_Meps(double m,double number); + /// Returns number farther from zero by smallest interval * 10^8 + double plus_2eps(double number); + /// Returns number closer to zero by smallest interval * 10^8 + double minus_2eps(double number); + + /// Returns false if the number is stored as NAN + bool not_nan(double number); + + // Matrix stuff + + /// Invert a matrix (Gauss-Jordan) + /** Generates error 303 L 9 for singular matrix + * No checking for memory limitations **/ + void invert(double **matrix, unsigned size, Error &error); + + /// Move a value from a fraction of one range to a fraction of another + /** am::rerange(0,1,0.3,0,100) will return 30. No checking to enforce + * that the value actually lies within the range is made. **/ + double rerange(double one_start, double one_end, double value, + double two_start, double two_end); +} + +#endif diff --git a/tools/pymol_asphere/src/spherical.cpp b/tools/pymol_asphere/src/spherical.cpp new file mode 100755 index 0000000000..f209975dd6 --- /dev/null +++ b/tools/pymol_asphere/src/spherical.cpp @@ -0,0 +1,145 @@ +/*************************************************************************** + spherical.cpp + W. Michael Brown + ------------------- + + Stuff for working spherical coordinates + + __________________________________________________________________________ + + Part of the Math Library + __________________________________________________________________________ + + begin : Tue Aug 29 2006 + copyright : (C) 2006 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#include "spherical.h" + +template class Ball; +template class Ball; +template ostream & operator<< (ostream &out, const Ball &t); +template ostream & operator<< (ostream &out, const Ball &t); +template istream & operator>> (istream &in, Ball &t); +template istream & operator>> (istream &in, Ball &t); + +// Empty construct. Not necessarily initialized to [0 0] +template +Ball::Ball() { +} + +// Assign theta and phi to the value +template +Ball::Ball(numtyp value) : theta(value), phi(value) { +} + +// Assignment Constructor +template +Ball::Ball(numtyp thet, numtyp ph) : theta(thet), phi(ph) { +} + +// Convert from cartesian +template +Ball::Ball(const ThreeD &pt) { + theta=atan2(pt.y,pt.x); + if (theta<0) + theta+=TWOPI; + phi=acos(pt.z/pt.norm()); +} + +template +numtyp & Ball::operator[](unsigned i) { + if (i==THETA) + return theta; + else + return phi; +} + +template +ostream & operator<< (ostream &out, const Ball &t) { + out << t.theta << " " << t.phi; + return out; +} + +template +istream & operator>> (istream &in, Ball &t) { + in >> t.theta >> t.phi; + return in; +} + +// Distance between two points (along arc) +template +numtyp Ball::dist(const Ball &two) const { + double dot=cPt(*this).dot(cPt(two)); + if (dot>1.0) + dot=1.0; + return acos(dot); +} + +// Distance squared between two points +template +numtyp Ball::dist2(const Ball &two) const { + numtyp d=dist(two); + return d*d; +} + +// Add both angles +template +void Ball::operator += (const Ball &two) { + theta+=two.theta; + phi+=two.phi; +} + +// Add to both angles +template +Ball Ball::operator + (const numtyp two) const { + return Ball(theta+two,phi+two); +} + +// Multiply both angles +template +Ball Ball::operator * (const numtyp two) const { + return Ball(theta*two,phi*two); +} + +// Divide both angles +template +void Ball::operator /= (const numtyp two) { + theta/=two; + phi/=two; +} + +// Add both angles +template +Ball Ball::operator + (const Ball &two) { + return Ball(theta+two.theta,phi+two.phi); +} + +// Move coordinates into array +template +void Ball::to_array(numtyp *array) { + *array=theta; + array++; + *array=phi; +} + +// Set coordinates from array +template +void Ball::from_array(numtyp *array) { + theta=*array; + array++; + phi=*array; +} + +// Returns 2 +template +unsigned Ball::dimensionality() { + return 2; +} + +// Returns true +template +bool Ball::check_bounds(numtyp min,numtyp max) { + return true; +} diff --git a/tools/pymol_asphere/src/spherical.h b/tools/pymol_asphere/src/spherical.h new file mode 100755 index 0000000000..f701625a79 --- /dev/null +++ b/tools/pymol_asphere/src/spherical.h @@ -0,0 +1,113 @@ +/*************************************************************************** + spherical.h + W. Michael Brown + ------------------- + + Stuff for working spherical coordinates + + __________________________________________________________________________ + + Part of the Math Library + __________________________________________________________________________ + + begin : Tue Aug 29 2006 + copyright : (C) 2006 by W. Michael Brown + email : wmbrown@sandia.gov + ***************************************************************************/ + +#ifndef SPHERICAL_H +#define SPHERICAL_H + +#include "miscm.h" +#include "m_constants.h" +#include "cartesian.h" +#include +#include +#include +#include +using namespace std; + +// Other coordinates +template class ThreeD; + +// Friends +template class Ball; +template +ostream & operator<< (ostream &out, const Ball &t); +template +istream & operator>> (istream &in, Ball &t); + +enum { THETA, ///<0 + PHI ///<1 +}; + +/// Two dimensional spherical coordinates on a unit sphere +/** The elements can be accessed directly .theta or .phi + * or by using the operator [] ( [THETA], [PHI] ) + * + * Input and output are overloaded for element I/O of the form "theta phi" + * <<, >> + **/ +template +class Ball { + public: + /// Empty construct. Not necessarily initialized to [0 0] + Ball(); + /// Assignment Constructor + Ball(numtyp theta, numtyp phi); + /// Assign theta and phi to the value + Ball(numtyp value); + /// Convert from cartesian + Ball(const ThreeD &pt); + + numtyp theta; + numtyp phi; + + numtyp &operator[](unsigned i); + + friend ostream & operator<< <>(ostream &out, const Ball &t); + friend istream & operator>> <>(istream &in, Ball &t); + + /// Add both angles + void operator += (const Ball &two); + /// Add to both angles + Ball operator + (const numtyp two) const; + /// Multiply both angles + Ball operator * (const numtyp two) const; + /// Divide both angles + void operator /= (const numtyp two); + /// Add both angles + Ball operator + (const Ball &two); + + /// Distance between two points (along arc) + /** \note The form of calculation used suffers from round off error + * when points are antipodal **/ + numtyp dist(const Ball &two) const; + /// Distance squared between two points (along arc) + /** \note The form of calculation used suffers from round off error + * when points are antipodal **/ + numtyp dist2(const Ball &two) const; + + /// Move coordinates into array + void to_array(numtyp *array); + /// Set coordinates from array + void from_array(numtyp *array); + + // -------------- Weird functions that help with coord templating + /// Returns 2 + unsigned dimensionality(); + // Returns true + bool check_bounds(numtyp min,numtyp max); + private: +}; + +///\var typedef Ball BallD +/// Double unit sphere +typedef Ball BallD; +///\var typedef Ball BallF +/// Float unit sphere +typedef Ball BallF; + + + +#endif