1 /*
2 * (c) 2004 Iowa State University
3 * see the LICENSE file in the top level directory
4 */
5
6 /*
7 Added RotateMoleculeGL routine and code to generate the angle strings and trackball - BMB July 2001
8 Corrected Create3DGLPICT for resolutions other than 72 dpi - BMB Feb 2002
9 */
10
11 #include <new>
12 #include <string.h>
13 #include <iostream>
14 #include <algorithm>
15 #include <wx/stdpaths.h>
16
17 #include "Globals.h"
18 #include "MoleculeData.h"
19 #include "Frame.h"
20 #include "SurfaceTypes.h"
21 #include "Math3D.h"
22 #include "InputData.h"
23 #include "MolDisplayWin.h"
24 #include "mpGLCanvas.h"
25 #include "mmp_gl.h"
26 #include "glf.h"
27 #include "Prefs.h"
28 #include "Progress.h"
29 #include "patterns.h"
30 #include "Files.h"
31 #include "Gradient.h"
32
33 // ----------------------------------------------------------------------------
34 // FUNCTION PROTOTYPES
35 // ----------------------------------------------------------------------------
36
37 void DrawPipeCylinder(float length, GLUquadric *quadric, unsigned int ncaps,
38 GLuint sphere_list, float radius, long quality);
39 void DrawPipeSpheres(const WinPrefs& Prefs, float length, float scale_factor,
40 bool is_masked, GLuint sphere_list);
41 void DashedQuadFromLine(const CPoint3D& pt1, const CPoint3D& pt2, float width,
42 float m[16], const CPoint3D& x_world, float offset,
43 GLuint length_anno_tex_id, const WinPrefs * Prefs,
44 bool draw_label = true);
45 void DrawAngleAnnotation(const CPoint3D *pt1, const CPoint3D *pt2,
46 const CPoint3D *pt3, const WinPrefs *Prefs,
47 GLuint length_anno_tex_id, CPoint3D *ambig_axis = NULL);
48 void CreateCylinderFromLine(GLUquadricObj * qobj, const CPoint3D & lineStart,
49 const CPoint3D & lineEnd, float lineWidth,
50 int nslices = 4, int nstacks = 1, bool cap = false);
51 void DrawRotationAxis(const CPoint3D & lineStart, const CPoint3D & lineEnd,
52 const int & order);
53 void DrawInversionPoint(void);
54 void DrawTranslucentPlane(const CPoint3D & origin, const CPoint3D & p1,
55 const CPoint3D & p2);
56 void DrawArrow(const float & length, const float & width, const int & quality);
57 void DrawSceneString(const float scale_factor, const float shift_x,
58 const float shift_y, const float shift_z,
59 const wxString& label);
60 void DrawString(const char *str);
61
62 // ----------------------------------------------------------------------------
63 // GLOBALS
64 // ----------------------------------------------------------------------------
65
66 const GLubyte stippleMask[128] = {
67 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
68 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00,
69 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
70 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00,
71 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
72 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00,
73 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
74 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00,
75 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
76 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00,
77 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
78 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00,
79 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
80 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00,
81 0xaa, 0xaa, 0xaa, 0xaa, 0x00, 0x00, 0x00, 0x00,
82 0x22, 0x22, 0x22, 0x22, 0x00, 0x00, 0x00, 0x00
83 };
84
85 GLfloat d_specular[] = {0.1f, 0.1f, 0.1f, 1.0};
86 GLfloat d_shininess = 1.0;
87 GLfloat d_diffuse[] = {0.02f,0.02f,0.02f,0.8f};
88 GLfloat d_ambient[] = {0.1f,0.1f,0.1f,0.8f};
89
90 GLfloat l_emissive[] = {0.0, 0.0, 0.0, 1.0};
91 GLfloat l_specular[] = {0.8f, 0.8f, 0.8f, 1.0};
92 GLfloat l_shininess = 80.0;
93 GLfloat l_diffuse[] = {0.2f,0.2f,0.2f,0.8f};
94 GLfloat l_ambient[] = {0.1f,0.1f,0.1f,0.8f};
95
96 // ----------------------------------------------------------------------------
97 // FUNCTION DEFINITIONS
98 // ----------------------------------------------------------------------------
99
100 /**
101 * This function frees any transparent geometries and any active display
102 * lists. It should be called when the view changes (which alters sort order
103 * of transparent surfaces), surfaces are updated, or the display list content
104 * has been altered. It does not remake any lists.
105 */
ReleaseLists()106 void MolDisplayWin::ReleaseLists() {
107 triangleCount = 0;
108
109 delete[] transpTriList;
110 transpTriList = NULL;
111 delete[] transpSortVertex;
112 transpSortVertex = NULL;
113 delete[] transpIndex;
114 transpIndex = NULL;
115
116 if (MainListActive) {
117 glDeleteLists(MainDisplayList, 1);
118 MainListActive = false;
119 }
120 if (SurfaceListActive) {
121 glDeleteLists(SurfaceDisplayList, 1);
122 SurfaceListActive = false;
123 }
124 }
125
126 /**
127 Display a transformation annotation, optionally consisting of Euler
128 angles and an outline of the trackball's footprint.
129 */
ShowRotation(bool ShowAngles,bool ShowTrackball)130 void MolDisplayWin::ShowRotation(bool ShowAngles, bool ShowTrackball) {
131
132 GLint matrixMode;
133 glGetIntegerv(GL_MATRIX_MODE, &matrixMode);
134 glMatrixMode(GL_PROJECTION);
135 glPushMatrix();
136 glLoadIdentity();
137 glMatrixMode(GL_MODELVIEW);
138 glPushMatrix();
139 glLoadIdentity();
140 wxRect DisplayRect = glCanvas->GetRect();
141 long hsize = DisplayRect.GetWidth();
142 long vsize = DisplayRect.GetHeight();
143 glScalef(2.0f / hsize, -2.0f / vsize, 1.0);
144 glTranslatef(-hsize / 2.0f, -vsize / 2.0f, 0.0);
145
146 glColor3fv(fg_color);
147
148 if (ShowAngles) {
149 char AngleString[50];
150 float psi, phi, theta;
151 MatrixToEulerAngles(MainData->TotalRotation, &psi, &phi, &theta);
152 sprintf((char *)AngleString, "%.2f, %.2f, %.2f, Scale:%.2f",
153 psi, phi, theta, MainData->WindowSize);
154
155 DrawStaticLabel(AngleString, 10, -20);
156 }
157
158 //Draw the trackball outline
159 if (ShowTrackball) {
160 wxPoint sphereCenter;
161 long sphereRadius;
162 sphereCenter.x = hsize / 2;
163 sphereCenter.y = vsize / 2;
164 if (sphereCenter.x >= sphereCenter.y)
165 sphereRadius = (long) (((float) sphereCenter.x) * 0.9);
166 else
167 sphereRadius = (long) (((float) sphereCenter.y) * 0.9);
168 long NumDivisions = (long) (20.0 * (1.0 + sphereRadius / 200.0));
169 float divarc = (2 * kPi) / NumDivisions;
170
171 glLineWidth(1);
172 glBegin(GL_LINE_LOOP);
173 glVertex3d(sphereCenter.x - sphereRadius, sphereCenter.y, 0.0);
174 for (int i=0; i<NumDivisions; i++) {
175 float x = sphereCenter.x - (sphereRadius*cos(i*divarc));
176 float y = sphereCenter.y + (sphereRadius*sin(i*divarc));
177 glVertex3d(x, y, 0.0);
178 }
179 glEnd();
180 }
181 glPopMatrix(); // GL_MODELVIEW
182 glMatrixMode(GL_PROJECTION);
183 glPopMatrix();
184 glMatrixMode(matrixMode);
185
186 }
187
188 /**
189 * Draws molecular structures, surfaces, and annotations.
190 * @pre Context, projection, and viewport have been set.
191 */
DrawGL(int do_shader)192 void MolDisplayWin::DrawGL(int do_shader) {
193 if (GLUpdateInProgress) return; //a previous call in the stack is already updating this window. Do not repeat.
194
195 GLUpdateInProgress = true; //set the flag to prevent calls from updating at the same time.
196
197 // Make (0, 0, WindowSize) the origin. Moves geometry away from camera.
198 glTranslatef(0, 0, -MainData->WindowSize);
199
200 if ((Prefs->ShowAtomicSymbolLabels() || Prefs->ShowAtomNumberLabels()) &&
201 (shader_program == 0 || do_shader)) {
202 glMatrixMode(GL_TEXTURE);
203 glPushMatrix();
204 glLoadIdentity();
205 glMatrixMode(GL_MODELVIEW);
206 DrawLabel();
207 glMatrixMode(GL_TEXTURE);
208 glPopMatrix();
209 glMatrixMode(GL_MODELVIEW);
210 }
211
212 glEnable(GL_LIGHTING);
213
214 // Rotate about the origin (the molecule's centroid).
215 glMultMatrixf((const GLfloat *) &(MainData->TotalRotation));
216
217 // Make the molecule's centroid the origin.
218 glTranslatef(-MainData->Centroid.x, -MainData->Centroid.y,
219 -MainData->Centroid.z);
220
221 // We're in a good space to draw the axis. We've rotated, but we haven't
222 // moved around for the molecule yet.
223 if (MainData->ShowAxis()) AddAxisGL();
224
225 if (GLEW_VERSION_2_0 && do_shader) {
226 GLint mode;
227 glGetIntegerv(GL_RENDER_MODE, &mode);
228 if (Prefs->GetShaderMode() && mode == GL_RENDER) {
229 glBindTexture(GL_TEXTURE_2D, depth_tex_id);
230 glUseProgram(shader_program);
231 }
232 }
233
234 //Draw the main molecular geometry
235 if (MainData->cFrame->NumAtoms > 0) {
236 if (MainListActive) {
237 glCallList(MainDisplayList);
238 } else { // build the main display list
239 // Suppress this temporarily because double- and triple-bond display
240 // requires a transformation depending on the current viewing
241 // transformation.
242 // MainDisplayList = glGenLists(1);
243 // glNewList(MainDisplayList, GL_COMPILE_AND_EXECUTE);
244 DrawMoleculeCoreGL();
245 // glEndList();
246 // MainListActive = true;
247 }
248 }
249
250 if (GLEW_VERSION_2_0) {
251 glUseProgram(0);
252 }
253
254 if (MainData->GetAnnotationCount() > 0) {
255 glMatrixMode(GL_TEXTURE);
256 glPushMatrix();
257 glLoadIdentity();
258 glMatrixMode(GL_MODELVIEW);
259 glLoadName(MMP_ANNOTATION);
260 glPushName(0);
261 std::vector<Annotation *>::const_iterator anno;
262 glDisable(GL_BLEND);
263 glDisable(GL_LIGHTING);
264 int anno_id = 0;
265 for (anno = MainData->Annotations.begin();
266 anno != MainData->Annotations.end(); anno++) {
267 glLoadName(anno_id + 1);
268 glColor3fv(fg_color);
269 (*anno)->draw(this);
270 anno_id++;
271 }
272 glPopName();
273 glEnable(GL_LIGHTING);
274 glMatrixMode(GL_TEXTURE);
275 glPopMatrix();
276 glMatrixMode(GL_MODELVIEW);
277 }
278
279 if (GLEW_VERSION_2_0 && Prefs->GetShaderMode() && do_shader) {
280 glBindTexture(GL_TEXTURE_2D, depth_tex_id);
281 glUseProgram(shader_program);
282 }
283
284 // Add any surfaces
285 Surface * lSurface = MainData->cFrame->SurfaceList;
286 haveTransparentSurfaces = false;
287
288 glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, l_emissive);
289 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, l_specular);
290 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 30.0f);
291
292 glEnable(GL_RESCALE_NORMAL);
293 int surf_id = 0;
294 while (lSurface) {
295 ++surf_id;
296 if (lSurface->GetVisibility()) {
297 if (!lSurface->isTransparent()) {
298 lSurface->Draw3DGL(MainData, Prefs, NULL);
299 } else {
300 haveTransparentSurfaces = true;
301 }
302 }
303 lSurface = lSurface->GetNextSurface();
304 }
305 glDisable(GL_RESCALE_NORMAL);
306
307 if (Prefs->ShowSymmetryOperators()) {
308 if (GLEW_VERSION_2_0) {
309 glUseProgram(0);
310 }
311 glMatrixMode(GL_TEXTURE);
312 glPushMatrix();
313 glLoadIdentity();
314 glMatrixMode(GL_MODELVIEW);
315 AddSymmetryOperators();
316 glMatrixMode(GL_TEXTURE);
317 glPopMatrix();
318 glMatrixMode(GL_MODELVIEW);
319 if (GLEW_VERSION_2_0 && Prefs->GetShaderMode() && do_shader) {
320 glBindTexture(GL_TEXTURE_2D, depth_tex_id);
321 glUseProgram(shader_program);
322 }
323 }
324
325 //Transparent surfaces have to be depth sorted and drawn last.
326 if (haveTransparentSurfaces) {
327 glEnable(GL_BLEND);
328 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
329 Surface * lSurface = MainData->cFrame->SurfaceList;
330 if (! transpTriList) {
331 long totalTriCount = 0;
332 while (lSurface) {
333 if (lSurface->GetVisibility() && lSurface->isTransparent()) {
334 totalTriCount += lSurface->getTriangleCount();
335 }
336 lSurface = lSurface->GetNextSurface();
337 }
338 transpTriList = new myGLTriangle[totalTriCount];
339 transpSortVertex = new float[totalTriCount];
340 transpIndex = new long[totalTriCount];
341 triangleCount = totalTriCount;
342
343 lSurface = MainData->cFrame->SurfaceList;
344 long triStartCount = 0;
345 while (lSurface) {
346 if (lSurface->GetVisibility() && lSurface->isTransparent()) {
347 triStartCount += lSurface->Draw3DGL(MainData, Prefs,
348 &(transpTriList[triStartCount]));
349 }
350 lSurface = lSurface->GetNextSurface();
351 }
352 for (int i=0; i<triangleCount; i++) transpIndex[i] = i;
353 SortTransparentTriangles();
354 }
355 DrawTransparentTriangles();
356 glDisable(GL_BLEND);
357 }
358
359 if (GLEW_VERSION_2_0) {
360 glUseProgram(0);
361 }
362
363 /* glPopMatrix(); // unrotated scene */
364
365 if (lasso_has_area) {
366 int canvas_width, canvas_height;
367
368 glEnable(GL_BLEND);
369 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
370
371 glMatrixMode(GL_MODELVIEW);
372 glPushMatrix();
373 glLoadIdentity();
374
375 glMatrixMode(GL_PROJECTION);
376 glPushMatrix();
377 glLoadIdentity();
378
379 glCanvas->GetClientSize(&canvas_width, &canvas_height);
380 gluOrtho2D(0, canvas_width, 0, canvas_height);
381
382 glColor4f(0.5f, 0.5f, 0.5f, 0.4f);
383 glRectf(lasso_start.x, lasso_start.y, lasso_end.x, lasso_end.y);
384
385 glDisable(GL_DEPTH_TEST);
386 glEnable(GL_LINE_STIPPLE);
387 glLineStipple(3, (GLushort) 43690);
388 glColor4f(fg_color[0], fg_color[1], fg_color[2], 0.9f);
389 glLineWidth(2.0f);
390 glBegin(GL_LINE_LOOP);
391 glVertex2f(lasso_start.x, lasso_start.y);
392 glVertex2f(lasso_end.x, lasso_start.y);
393 glVertex2f(lasso_end.x, lasso_end.y);
394 glVertex2f(lasso_start.x, lasso_end.y);
395 glEnd();
396 glLineWidth(1.0f);
397 glDisable(GL_LINE_STIPPLE);
398 glEnable(GL_DEPTH_TEST);
399
400 glPopMatrix();
401
402 glMatrixMode(GL_MODELVIEW);
403 glPopMatrix();
404
405 glDisable(GL_BLEND);
406 }
407
408 glDisable(GL_LIGHTING);
409
410 if (do_rotate_annotation) {
411 glMatrixMode(GL_TEXTURE);
412 glPushMatrix();
413 glLoadIdentity();
414 glMatrixMode(GL_MODELVIEW);
415 ShowRotation(Prefs->GetShowAngles() && !rotate_timer.IsRunning(),
416 !rotate_timer.IsRunning());
417 glMatrixMode(GL_TEXTURE);
418 glPopMatrix();
419 glMatrixMode(GL_MODELVIEW);
420 }
421 GLUpdateInProgress = false; //reset the flag to allow updates.
422 }
423
draw(const MolDisplayWin * win) const424 void AnnotationLength::draw(const MolDisplayWin * win) const {
425 MoleculeData * maindata = win->GetData();
426 Frame * cFrame = maindata->GetCurrentFramePtr();
427 WinPrefs * Prefs = win->GetPrefs();
428 //validate the atom references for this frame
429 if ((atoms[0] < 0)||(atoms[0] >= cFrame->GetNumAtoms())) return;
430 if ((atoms[1] < 0)||(atoms[1] >= cFrame->GetNumAtoms())) return;
431
432 GLdouble modelview[16];
433 GLdouble proj[16];
434 GLint viewport[4];
435 glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
436 glGetDoublev(GL_PROJECTION_MATRIX, proj);
437 glGetIntegerv(GL_VIEWPORT, viewport);
438
439 CPoint3D lookat_eye = CPoint3D(0.0f, 0.0f, 1.0f);
440 CPoint3D up_eye = CPoint3D(0.0f, 1.0f, 0.0f);
441 CPoint3D lookat_world;
442 CPoint3D up_world;
443 CPoint3D r;
444 float m[16];
445 Matrix4D mv_inv;
446
447 // What we want to do here is make the annotation always face the
448 // viewer. The computations below compute a rotation matrix to align
449 // to the camera. This only needs to be computed once for all
450 // annotations since the camera doesn't change positions.
451
452 // Invert just the rotation portion of the modelview matrix. This is
453 // much faster than inverting an arbitrary matrix.
454 mv_inv[0][0] = modelview[0];
455 mv_inv[0][1] = modelview[4];
456 mv_inv[0][2] = modelview[8];
457 mv_inv[0][3] = 0.0f;
458 mv_inv[1][0] = modelview[1];
459 mv_inv[1][1] = modelview[5];
460 mv_inv[1][2] = modelview[9];
461 mv_inv[1][3] = 0.0f;
462 mv_inv[2][0] = modelview[2];
463 mv_inv[2][1] = modelview[6];
464 mv_inv[2][2] = modelview[10];
465 mv_inv[2][3] = 0.0f;
466 mv_inv[3][0] = 0.0f;
467 mv_inv[3][1] = 0.0f;
468 mv_inv[3][2] = 0.0f;
469 mv_inv[3][3] = 1.0f;
470
471 // Transform the eye space vectors to world coordinates, and find
472 // a third vector to form a basis set.
473 Rotate3DPt(mv_inv, lookat_eye, &lookat_world);
474 Rotate3DPt(mv_inv, up_eye, &up_world);
475 CrossProduct3D(&lookat_world, &up_world, &r);
476
477 m[0] = r.x;
478 m[1] = r.y;
479 m[2] = r.z;
480 m[3] = 0.0f;
481
482 m[4] = up_world.x;
483 m[5] = up_world.y;
484 m[6] = up_world.z;
485 m[7] = 0.0f;
486
487 m[8] = lookat_world.x;
488 m[9] = lookat_world.y;
489 m[10] = lookat_world.z;
490 m[11] = 0.0f;
491
492 m[12] = m[13] = m[14] = 0.0f;
493 m[15] = 1.0f;
494
495 // x_world will indicate what vector in world coordinates will effect
496 // a direction in the eye's x direction. This is the direction in
497 // which the length label will appear.
498 CPoint3D x_eye;
499 CPoint3D x_world;
500 x_eye.x = 1.0f;
501 x_eye.y = 0.0f;
502 x_eye.z = 0.0f;
503 Rotate3DPt(mv_inv, x_eye, &x_world);
504
505 int bond_id;
506 float bond_size;
507
508 bond_id = cFrame->BondExists(atoms[0], atoms[1]);
509
510 // If a bond exists between the two atoms, we need to push out the
511 // length label accordingly.
512 GLdouble BondSize = Prefs->GetQD3DBondWidth();
513 if (bond_id > -1) {
514 bond_size = BondSize / MAX(cFrame->GetBondOrder(bond_id), 1.0f) *
515 3.5f * cFrame->GetBondOrder(bond_id) / 2.0f;
516 if (cFrame->GetBondOrder(bond_id) > 1) {
517 bond_size *= 1.5;
518 }
519 } else {
520 bond_size = 0.0f;
521 }
522
523 CPoint3D pt1, pt2;
524 cFrame->GetAtomPosition(atoms[0], pt1);
525 cFrame->GetAtomPosition(atoms[1], pt2);
526
527 // Draw the dashed line and label.
528 DashedQuadFromLine(pt1, pt2, BondSize * 0.25, m, x_world, bond_size,
529 win->GetLengthTexId(), Prefs);
530 }
531
draw(const MolDisplayWin * win) const532 void AnnotationAngle::draw(const MolDisplayWin * win) const {
533
534 MoleculeData *maindata = win->GetData();
535 Frame *cFrame = maindata->GetCurrentFramePtr();
536
537 // Validate the atom references for this frame
538 if ((atoms[0] < 0)||(atoms[0] >= cFrame->GetNumAtoms())) return;
539 if ((atoms[1] < 0)||(atoms[1] >= cFrame->GetNumAtoms())) return;
540 if ((atoms[2] < 0)||(atoms[2] >= cFrame->GetNumAtoms())) return;
541
542 CPoint3D atom1_pos, atom2_pos, atom3_pos;
543
544 cFrame->GetAtomPosition(atoms[0], atom1_pos);
545 cFrame->GetAtomPosition(atoms[1], atom2_pos);
546 cFrame->GetAtomPosition(atoms[2], atom3_pos);
547
548 DrawAngleAnnotation(&atom1_pos, &atom2_pos, &atom3_pos, win->GetPrefs(),
549 win->GetLengthTexId());
550
551 }
552
draw(const MolDisplayWin * win) const553 void AnnotationDihedral::draw(const MolDisplayWin * win) const {
554
555 MoleculeData *maindata = win->GetData();
556 Frame *cFrame = maindata->GetCurrentFramePtr();
557
558 //validate the atom references for this frame
559 if ((atoms[0] < 0)||(atoms[0] >= cFrame->GetNumAtoms())) return;
560 if ((atoms[1] < 0)||(atoms[1] >= cFrame->GetNumAtoms())) return;
561 if ((atoms[2] < 0)||(atoms[2] >= cFrame->GetNumAtoms())) return;
562 if ((atoms[3] < 0)||(atoms[3] >= cFrame->GetNumAtoms())) return;
563
564 CPoint3D atom1_pos;
565 CPoint3D atom2_pos;
566 CPoint3D atom3_pos;
567 CPoint3D atom4_pos;
568 CPoint3D vec[2];
569 CPoint3D vec2;
570 CPoint3D normal1;
571 CPoint3D normal2;
572 CPoint3D binormal1;
573 CPoint3D binormal2;
574 Matrix4D plane2xy;
575 float angle;
576 GLfloat color[4];
577
578 glGetFloatv(GL_CURRENT_COLOR, color);
579
580 cFrame->GetAtomPosition(atoms[0], atom1_pos);
581 cFrame->GetAtomPosition(atoms[1], atom2_pos);
582 cFrame->GetAtomPosition(atoms[2], atom3_pos);
583 cFrame->GetAtomPosition(atoms[3], atom4_pos);
584
585 // The first three atoms (1, 2, 3) form one plane. The last three
586 // (2, 3, 4) form the second plane. We want to find the angle
587 // between those two planes, for that is the dihedral angle.
588
589 // Find the first plane's normal by finding the vectors between its
590 // points and calculating the normal from them.
591 vec[0] = atom1_pos - atom3_pos;
592 Normalize3D(&vec[0]);
593
594 vec2 = atom2_pos - atom3_pos;
595 Normalize3D(&vec2);
596
597 CrossProduct3D(&vec[0], &vec2, &normal1);
598
599 // Find the second plane. The two planes share a vector between
600 // atoms 2 and 3, so we just reuse that from the first plane.
601 vec[1] = atom4_pos - atom3_pos;
602 Normalize3D(&vec[1]);
603
604 CrossProduct3D(&vec[1], &vec2, &normal2);
605
606 // Okay, we've drawn the planes. Now we want to draw an angle
607 // annotation between them. The angle annotation should span
608 // from the midpoint on one half-circle's arc to the midpoint on
609 // the other half circles arc. These can be found by considering
610 // the vector between atoms 2 and 3 as the tangent vector, and
611 // find the the binormal as the crossproduct of tangent and normal.
612 CrossProduct3D(&vec2, &normal1, &binormal1);
613 CrossProduct3D(&vec2, &normal2, &binormal2);
614
615 Normalize3D(&binormal1);
616 Normalize3D(&binormal2);
617
618 // Figure out the actual points from the vector displacements.
619 CPoint3D pt1 = atom3_pos + binormal1;
620 CPoint3D pt3 = atom3_pos + binormal2;
621
622 Matrix4D m;
623 CPoint3D pt1_eye;
624 CPoint3D pt3_eye;
625 int id;
626
627 glGetFloatv(GL_MODELVIEW_MATRIX, (float *) m);
628 Rotate3DPt(m, pt1, &pt1_eye);
629 Rotate3DPt(m, pt3, &pt3_eye);
630
631 id = pt3_eye.z < pt1_eye.z;
632
633 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
634 glEnable(GL_BLEND);
635 glDepthMask(GL_FALSE);
636 for (int j = 0; j <= 1; j++, id = (id + 1) % 2) {
637 // Make this plane look like the x-y plane for easier circle
638 // drawing. The vector from atom3 to atom2 should look like
639 // the positive x-axis.
640 SetPlaneRotation(plane2xy, vec2, vec[id]);
641
642 glPushMatrix();
643 glTranslatef(atom3_pos.x, atom3_pos.y, atom3_pos.z);
644 glMultMatrixf((GLfloat *) plane2xy);
645
646 glColor4f(color[0] - 0.3f, color[1] - 0.3f, color[2] - 0.3f, 0.2f);
647 glBegin(GL_TRIANGLE_FAN);
648 glVertex3f(0.0f, 0.0f, 0.0f);
649 for (angle = 0.0f; angle <= 3.1416f; angle += 0.01f) {
650 glVertex3f(cos(angle), sin(angle), 0.0f);
651 }
652 glVertex3f(0.0f, 0.0f, 0.0f);
653 glEnd();
654
655 glColor4f(color[0] - 0.6f, color[1] - 0.06f, color[2] - 0.06f, 0.2f);
656 glLineWidth(2.0f);
657 glBegin(GL_LINES);
658 glVertex3f(-1.0f, 0.0f, 0.0f);
659 glVertex3f(1.0f, 0.0f, 0.0f);
660 glEnd();
661 glLineWidth(1.0f);
662
663 glPopMatrix();
664 }
665
666 glDepthMask(GL_TRUE);
667 glDisable(GL_BLEND);
668
669 glColor4f(color[0], color[1], color[2], 1.0f);
670 DrawAngleAnnotation(&pt1, &atom3_pos, &pt3, win->GetPrefs(),
671 win->GetLengthTexId(), &normal1);
672
673 }
674
draw(const MolDisplayWin * win) const675 void AnnotationMarker::draw(const MolDisplayWin * win) const {
676
677 float m[16];
678 float angle;
679 float ca1, sa1;
680 float ca2, sa2;
681 float invert1[3];
682 float invert2[3];
683 CPoint3D atom_pos;
684 float radius;
685
686 MoleculeData *maindata = win->GetData();
687 Frame *cFrame = maindata->GetCurrentFramePtr();
688 WinPrefs * Prefs = win->GetPrefs();
689
690 // Validate the atom reference for this frame.
691 if (atoms[0] < 0 || atoms[0] >= cFrame->GetNumAtoms()) return;
692
693 float AtomScale = Prefs->GetAtomScale();
694 long curAtomType = cFrame->GetAtomType(atoms[0]) - 1;
695
696 radius = AtomScale * Prefs->GetAtomSize(curAtomType);
697
698 cFrame->GetAtomPosition(atoms[0], atom_pos);
699
700 glPushMatrix();
701 glTranslatef(atom_pos.x, atom_pos.y, atom_pos.z);
702
703 // RGBColor * BackgroundColor = Prefs->GetBackgroundColorLoc();
704
705 // invert1[0] = 1.0f - BackgroundColor->red / 65536.0f;
706 // invert1[1] = 1.0f - BackgroundColor->green / 65536.0f;
707 // invert1[2] = 1.0f - BackgroundColor->blue / 65536.0f;
708 invert1[0] = invert1[1] = invert1[2] = 0.0f;
709
710 glGetFloatv(GL_MODELVIEW_MATRIX, m);
711 glLineWidth(3.0f);
712
713 invert2[0] = 0.0f;
714 invert2[1] = 1.0f;
715 invert2[2] = 0.0f;
716
717 glBegin(GL_QUAD_STRIP);
718 for (angle = 0.0f; angle <= 6.24f; angle += 0.1f) {
719 ca1 = cos(angle) * (radius + 0.03f);
720 sa1 = sin(angle) * (radius + 0.03f);
721 ca2 = cos(angle) * (radius + 0.08f);
722 sa2 = sin(angle) * (radius + 0.08f);
723 glColor3fv(invert1);
724 glVertex3f(ca1 * m[0] + sa1 * m[1],
725 ca1 * m[4] + sa1 * m[5],
726 ca1 * m[8] + sa1 * m[9]);
727 glColor3fv(invert2);
728 glVertex3f(ca2 * m[0] + sa2 * m[1],
729 ca2 * m[4] + sa2 * m[5],
730 ca2 * m[8] + sa2 * m[9]);
731 }
732 ca1 = cos(0.0f) * (radius + 0.03f);
733 sa1 = sin(0.0f) * (radius + 0.03f);
734 ca2 = cos(0.0f) * (radius + 0.08f);
735 sa2 = sin(0.0f) * (radius + 0.08f);
736 glColor3fv(invert1);
737 glVertex3f(ca1 * m[0] + sa1 * m[1],
738 ca1 * m[4] + sa1 * m[5],
739 ca1 * m[8] + sa1 * m[9]);
740 glColor3fv(invert2);
741 glVertex3f(ca2 * m[0] + sa2 * m[1],
742 ca2 * m[4] + sa2 * m[5],
743 ca2 * m[8] + sa2 * m[9]);
744 glEnd();
745
746 glBegin(GL_QUAD_STRIP);
747 for (angle = 0.0f; angle <= 6.24f; angle += 0.1f) {
748 ca1 = cos(angle) * (radius + 0.08f);
749 sa1 = sin(angle) * (radius + 0.08f);
750 ca2 = cos(angle) * (radius + 0.13f);
751 sa2 = sin(angle) * (radius + 0.13f);
752 glColor3fv(invert2);
753 glVertex3f(ca1 * m[0] + sa1 * m[1],
754 ca1 * m[4] + sa1 * m[5],
755 ca1 * m[8] + sa1 * m[9]);
756 glColor3fv(invert1);
757 glVertex3f(ca2 * m[0] + sa2 * m[1],
758 ca2 * m[4] + sa2 * m[5],
759 ca2 * m[8] + sa2 * m[9]);
760 }
761 ca1 = cos(0.0f) * (radius + 0.08f);
762 sa1 = sin(0.0f) * (radius + 0.08f);
763 ca2 = cos(0.0f) * (radius + 0.13f);
764 sa2 = sin(0.0f) * (radius + 0.13f);
765 glColor3fv(invert2);
766 glVertex3f(ca1 * m[0] + sa1 * m[1],
767 ca1 * m[4] + sa1 * m[5],
768 ca1 * m[8] + sa1 * m[9]);
769 glColor3fv(invert1);
770 glVertex3f(ca2 * m[0] + sa2 * m[1],
771 ca2 * m[4] + sa2 * m[5],
772 ca2 * m[8] + sa2 * m[9]);
773 glEnd();
774 glLineWidth(1.0f);
775 glPopMatrix();
776 }
777
DrawStaticLabel(const char * label,GLfloat x,GLfloat y)778 void MolDisplayWin::DrawStaticLabel(const char *label, GLfloat x, GLfloat y) {
779 int canvasWidth, canvasHeight;
780
781 glCanvas->GetSize(&canvasWidth, &canvasHeight);
782 glColor3fv(fg_color);
783
784 GLint matrixMode;
785 glGetIntegerv(GL_MATRIX_MODE, &matrixMode);
786
787 glMatrixMode(GL_PROJECTION);
788 glPushMatrix();
789 glLoadIdentity();
790 gluOrtho2D(0, canvasWidth, 0, canvasHeight);
791
792 glMatrixMode(GL_MODELVIEW);
793 glPushMatrix();
794 glLoadIdentity();
795
796 if (x < 0) x = canvasWidth + x;
797 if (y < 0) y = canvasHeight + y;
798
799 glTranslatef(x, y, 0.0f);
800 glScalef(300, 300, 1);
801
802 glLoadName(MMP_NULL);
803
804 // glDisable(GL_DEPTH_TEST);
805 glfStartBitmapDrawing();
806
807 glBlendFunc(GL_SRC_ALPHA, GL_DST_ALPHA);
808 glEnable(GL_BLEND);
809 glfDrawBMaskString(label);
810 glDisable(GL_BLEND);
811
812 glfStopBitmapDrawing();
813 // glEnable(GL_DEPTH_TEST);
814
815 glFlush();
816
817 glPopMatrix(); // GL_MODELVIEW
818
819 glMatrixMode(GL_PROJECTION);
820 glPopMatrix();
821
822 glMatrixMode(matrixMode);
823 }
824
DrawLabel()825 void MolDisplayWin::DrawLabel() {
826 Frame * lFrame=MainData->cFrame;
827 mpAtom * lAtoms = lFrame->Atoms;
828 long NumAtoms = lFrame->NumAtoms;
829 float AtomScale = Prefs->GetAtomScale();
830 GLdouble BondSize = Prefs->GetQD3DBondWidth();
831 float LabelSize = Prefs->GetAtomLabelSize();
832 glfStringCentering(true);
833
834 wxString atomLabel;
835 long CurrentAtomType;
836 CPoint3D origPt, transPt;
837 if (!Prefs->DrawWireFrame() || Prefs->ColorBondHalves()) {
838 glLoadName(MMP_ATOM);
839 glPushName(0);
840 for (long iatom=0; iatom<NumAtoms; iatom++) {
841 if (lAtoms[iatom].GetInvisibility()) continue; //Atom is invisible so skip
842
843 atomLabel.Clear();
844
845 CurrentAtomType = lAtoms[iatom].GetType() - 1;
846
847 //!!! retrieve atom label
848 if ( Prefs->ShowAtomicSymbolLabels() )
849 Prefs->GetAtomLabel(CurrentAtomType, atomLabel);
850
851 if (Prefs->ShowAtomNumberLabels() ) {
852 wxString tmpStr;
853
854 tmpStr.Printf(wxT("%ld"), iatom+1);
855 atomLabel.Append(tmpStr);
856 }
857
858 float radius;
859 if (!Prefs->DrawWireFrame()) radius = AtomScale*Prefs->GetAtomSize(CurrentAtomType);
860 else radius = BondSize;
861
862 if (radius<0.01) continue; //skip really small spheres
863
864 RGBColor * AtomColor = Prefs->GetAtomColorLoc(CurrentAtomType);
865 float red, green, blue;
866 red = AtomColor->red/65536.0;
867 green = AtomColor->green/65536.0;
868 blue = AtomColor->blue/65536.0;
869
870 // Find atom's position in world coordinates. Move there for
871 // drawing label screen-aligned.
872 origPt = lAtoms[iatom].Position;
873 origPt -= MainData->Centroid;
874 Rotate3DPt(MainData->TotalRotation, origPt, &transPt);
875
876 glPushMatrix();
877 glTranslatef(transPt.x, transPt.y, transPt.z+(radius+0.01));
878
879 glColor3f(1-red, 1-green, 1-blue);
880 glScalef((0.1+0.08*radius)*LabelSize, (0.1+0.08*radius)*LabelSize, 1);
881 glLoadName(iatom+1);
882 glfDrawSolidString(atomLabel.mb_str(wxConvUTF8));
883 glPopMatrix();
884 }
885 glPopName();
886 }
887 glfStringCentering(false);
888 }
889
SortTransparentTriangles(void)890 void MolDisplayWin::SortTransparentTriangles(void) {
891 CPoint3D tempV;
892 for (int i=0; i<triangleCount; i++) {
893 Rotate3DPt(MainData->TotalRotation, transpTriList[i].v1,
894 &(tempV));
895 transpSortVertex[i] = tempV.z; //We only need to save the z coordinate to sort
896 }
897 bool done = false;
898 long maxcount = triangleCount-1;
899 while (!done) {
900 done = true;
901 for (int i=0; i<maxcount; i++) {
902 long tempi = transpIndex[i];
903 long tempj = transpIndex[i+1];
904 if (transpSortVertex[tempi] > transpSortVertex[tempj]) {
905 transpIndex[i] = tempj;
906 transpIndex[i+1] = tempi;
907 done = false;
908 }
909 }
910 maxcount--;
911 }
912 }
913
DrawTransparentTriangles(void)914 void MolDisplayWin::DrawTransparentTriangles(void) {
915 glBegin(GL_TRIANGLES);
916 for (long itri=0; itri<triangleCount; itri++) {
917 long itriIndex = transpIndex[itri];
918
919 glColor4f(transpTriList[itriIndex].r1, transpTriList[itriIndex].g1,
920 transpTriList[itriIndex].b1, transpTriList[itriIndex].a1);
921 glNormal3f(transpTriList[itriIndex].n1.x,
922 transpTriList[itriIndex].n1.y, transpTriList[itriIndex].n1.z);
923 glVertex3d(transpTriList[itriIndex].v1.x, transpTriList[itriIndex].v1.y,
924 transpTriList[itriIndex].v1.z);
925
926 glColor4f(transpTriList[itriIndex].r2, transpTriList[itriIndex].g2,
927 transpTriList[itriIndex].b2, transpTriList[itriIndex].a2);
928 glNormal3f(transpTriList[itriIndex].n2.x,
929 transpTriList[itriIndex].n2.y, transpTriList[itriIndex].n2.z);
930 glVertex3d(transpTriList[itriIndex].v2.x, transpTriList[itriIndex].v2.y,
931 transpTriList[itriIndex].v2.z);
932
933 glColor4f(transpTriList[itriIndex].r3, transpTriList[itriIndex].g3,
934 transpTriList[itriIndex].b3, transpTriList[itriIndex].a3);
935 glNormal3f(transpTriList[itriIndex].n3.x,
936 transpTriList[itriIndex].n3.y, transpTriList[itriIndex].n3.z);
937 glVertex3d(transpTriList[itriIndex].v3.x, transpTriList[itriIndex].v3.y,
938 transpTriList[itriIndex].v3.z);
939 }
940 glEnd(); //End of triangle creation
941 }
942
DrawMoleculeCoreGL(void)943 void MolDisplayWin::DrawMoleculeCoreGL(void) {
944
945 GLUquadric *core_obj;
946 core_obj = gluNewQuadric();
947
948 // gluQuadricDrawStyle(qobj, GLU_FILL); //or GLU_LINE
949 /* gluQuadricOrientation(core_obj, GLU_OUTSIDE); */
950 /* gluQuadricNormals(core_obj, GLU_SMOOTH); //GLU_FLAT GLU_NONE */
951
952 Frame *lFrame = MainData->cFrame;
953 mpAtom *lAtoms = lFrame->Atoms;
954 Bond *lBonds = lFrame->Bonds;
955 long NumAtoms = lFrame->NumAtoms;
956 long NumBonds = lFrame->NumBonds;
957 float AtomScale = Prefs->GetAtomScale();
958 float Quality = Prefs->GetQD3DAtomQuality();
959 // For really large systems the size of the individual atoms and bonds gets small so there doesn't
960 // seem to be a need for the quality to remain high. This is an attempt to silently degrade the display quality
961 // to improve performance.
962 // Could disable the reduction when printing if this ends up causing a noticable issue
963 int QualReduce = MainData->WindowSize;
964 QualReduce /= 15; //truncate the division
965 Quality -= QualReduce;
966 Quality = MAX(Quality, 2); //Don't go below 2
967
968 bool draw_subdued;
969 glEnable(GL_RESCALE_NORMAL);
970 if (!Prefs->DrawWireFrame()) {
971 glLoadName(MMP_ATOM);
972 glPushName(0);
973 for (long iatom=0; iatom<NumAtoms; iatom++) {
974
975 if (lAtoms[iatom].GetInvisibility()) continue; //Atom is invisible so skip
976 if (Prefs->ShowEFPWireFrame()&&lAtoms[iatom].IsEffectiveFragment()) continue;
977 long curAtomType = lAtoms[iatom].GetType() - 1;
978
979 float radius = AtomScale * Prefs->GetAtomSize(curAtomType);
980 if (radius < 0.01) continue; //skip really small spheres
981
982 /* glPushMatrix(); // atom center */
983 /* glCallList(sphere_list); */
984 /* glPopMatrix(); // atom center */
985
986 glPushMatrix(); // molecule origin
987 glTranslatef(lAtoms[iatom].Position.x,
988 lAtoms[iatom].Position.y,
989 lAtoms[iatom].Position.z);
990
991 if (InSymmetryEditMode() && !lAtoms[iatom].IsSymmetryUnique()) {
992 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, d_specular);
993 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, d_shininess);
994 glLoadName(0);
995 glColor3fv(fg_color);
996 glEnable(GL_POLYGON_STIPPLE);
997 glPolygonStipple(stippleMask);
998 gluSphere(core_obj, radius*1.01, (long)(1.5*Quality), (long)(Quality));
999 glDisable(GL_POLYGON_STIPPLE);
1000 glPopMatrix(); // molecule origin
1001 continue;
1002 }
1003
1004 glLoadName(iatom + 1);
1005 draw_subdued = mHighliteState && !lAtoms[iatom].GetSelectState();
1006
1007 if (draw_subdued) {
1008 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, d_specular);
1009 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, d_shininess);
1010 } else {
1011 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, l_specular);
1012 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, l_shininess);
1013 }
1014
1015 if (Prefs->Show2DPattern()) {
1016 short patternindex = Prefs->GetAtomPattern(lAtoms[iatom].GetType()-1);
1017 //The 0th pattern is assumed to be solid so no need to draw
1018 if ((patternindex>0)&&(patternindex<kNumPatterns)) {
1019 glColor3f(0.0f, 0.0f, 0.0f);
1020 glEnable(GL_POLYGON_STIPPLE);
1021 glPolygonStipple(atomMaskPatterns[patternindex]);
1022 glPushMatrix(); // atom center
1023 glScalef(radius, radius, radius);
1024 glCallList(sphere_list);
1025 glPopMatrix(); // atom center
1026 glDisable(GL_POLYGON_STIPPLE);
1027 }
1028 }
1029
1030 if (draw_subdued) {
1031 glColor3f(0.0f,0.0f,0.0f);
1032 glEnable(GL_POLYGON_STIPPLE);
1033 glPolygonStipple(stippleMask);
1034 glPushMatrix(); // atom center
1035 glScalef(radius * 1.01f, radius * 1.01f, radius * 1.01f);
1036 glCallList(sphere_list);
1037 glPopMatrix(); // atom center
1038 glDisable(GL_POLYGON_STIPPLE);
1039 }
1040
1041 Prefs->ChangeColorAtomColor(curAtomType+1);
1042
1043 if (draw_subdued) {
1044 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, d_specular);
1045 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, d_shininess);
1046 } else {
1047 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, l_specular);
1048 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, l_shininess);
1049 }
1050
1051 glPushMatrix(); // atom center
1052 glScalef(radius, radius, radius);
1053 glCallList(sphere_list);
1054 glPopMatrix(); // atom center
1055
1056 if (InEditMode() && show_bond_sites) {
1057 glPushMatrix();
1058 DrawBondingSites(iatom, radius, core_obj);
1059 glPopMatrix();
1060 }
1061
1062 glPopMatrix(); // molecular origin
1063 }
1064 glPopName();
1065 }
1066
1067 GLdouble modelview[16];
1068 GLdouble proj[16];
1069 GLint viewport[4];
1070
1071 glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
1072 glGetDoublev(GL_PROJECTION_MATRIX, proj);
1073 glGetIntegerv(GL_VIEWPORT, viewport);
1074
1075 // bonds as cylinders
1076 // In wireframe mode with bonds colored by atom color we simply scink the
1077 // atom radius to the bond size and get a nice rounded end cap. If bonds
1078 // are not colored by atom color then the sphere is skipped and a simple
1079 // disk closes off the cylinder
1080 glLoadName(MMP_BOND);
1081 glPushName(0);
1082 for (long ibond = 0; ibond < NumBonds; ibond++) {
1083
1084 if (!InSymmetryEditMode() ||
1085 (lAtoms[lBonds[ibond].Atom1].IsSymmetryUnique() &&
1086 lAtoms[lBonds[ibond].Atom2].IsSymmetryUnique())) {
1087 glLoadName(ibond + 1);
1088 } else {
1089 glLoadName(0);
1090 }
1091
1092 if (mHighliteState && !lBonds[ibond].GetSelectState()) {
1093 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, d_specular);
1094 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, d_shininess);
1095 } else {
1096 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, l_specular);
1097 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, l_shininess);
1098 }
1099
1100 DrawBond(lBonds[ibond], lAtoms[lBonds[ibond].Atom1],
1101 lAtoms[lBonds[ibond].Atom2], *Prefs, core_obj,
1102 modelview, proj, viewport, sphere_list,
1103 mHighliteState, InSymmetryEditMode());
1104
1105 }
1106 glPopName();
1107 glDisable(GL_RESCALE_NORMAL);
1108
1109 glLoadName(0); //only atoms and bonds are selectable
1110 //so give a NULL name value to the rest of the geometries
1111
1112 if (MainData->GetDrawMode() && lFrame->Vibs) { //Add the current normal mode, if active
1113 float VectorScale = Prefs->GetVectorScale();
1114
1115 long cmode = (lFrame->NumAtoms)*(lFrame->Vibs->CurrentMode);
1116
1117 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, l_specular);
1118 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 30.0f);
1119 Prefs->ChangeColorVectorColor();
1120
1121 CPoint3D NormStart = CPoint3D(0.0f, 0.0f, 1.0f);
1122 for (long iatom=0; iatom<NumAtoms; iatom++) {
1123
1124 if (lAtoms[iatom].GetInvisibility()) continue; //Atom is invisible so skip
1125 CPoint3D NMode = (lFrame->Vibs->NormMode[iatom + cmode]);
1126 NMode *= VectorScale;
1127 float length = NMode.Magnitude();
1128 if (length > 0.1) {
1129 CPoint3D NModeVector;
1130 Matrix4D rotMat;
1131 NModeVector.x = NMode.x/length;
1132 NModeVector.y = NMode.y/length;
1133 NModeVector.z = NMode.z/length;
1134 float VectorWidth = (0.03 + 0.005*length);
1135 //Set up vectors for the shaft of the normal mode
1136 CPoint3D VStart;
1137 long curAtomType = lAtoms[iatom].Type - 1;
1138
1139 float radius;
1140 if (Prefs->DrawWireFrame()) radius = 0.0;
1141 else radius = AtomScale*Prefs->GetAtomSize(curAtomType);
1142 VStart.x = lAtoms[iatom].Position.x + NModeVector.x*0.95*radius;
1143 VStart.y = lAtoms[iatom].Position.y + NModeVector.y*0.95*radius;
1144 VStart.z = lAtoms[iatom].Position.z + NModeVector.z*0.95*radius;
1145 float HeadRadius = 2 * VectorWidth;
1146 if (2*HeadRadius > length) HeadRadius = length/2.0;
1147 // float HeadRatio = (length-HeadRadius)/length;
1148 GLfloat ShaftLength = length - HeadRadius;
1149
1150 SetRotationMatrix(rotMat, &NormStart, &NModeVector);
1151 rotMat[3][0] = VStart.x;
1152 rotMat[3][1] = VStart.y;
1153 rotMat[3][2] = VStart.z;
1154 glPushMatrix();
1155 glMultMatrixf((const GLfloat *) &rotMat);
1156
1157 gluCylinder(core_obj, VectorWidth, VectorWidth, ShaftLength, (long)(Quality), (long)(0.5*Quality));
1158 glPopMatrix();
1159 rotMat[3][0] = VStart.x + NModeVector.x * ShaftLength;
1160 rotMat[3][1] = VStart.y + NModeVector.y * ShaftLength;
1161 rotMat[3][2] = VStart.z + NModeVector.z * ShaftLength;
1162 glPushMatrix();
1163 glMultMatrixf((const GLfloat *) &rotMat);
1164 gluDisk(core_obj, 0.0, 2*VectorWidth, (long)(Quality), 2);
1165 gluCylinder(core_obj, 2*VectorWidth, 0.0, HeadRadius, (long)(Quality), 3);
1166 glPopMatrix();
1167 }
1168 }
1169 }
1170
1171 if (MainData->GradientVectorAvailable() && Prefs->DisplayGradient()) { //Add the gradient arrows, if active
1172 float VectorScale = Prefs->GetGradientScale();
1173 // float GradMax = sqrt(1.0/lFrame->Gradient->GetMaximum());
1174 // float GradMax = lFrame->Gradient->GetMaximum()*10000.0;
1175 float GradMax = 0.01/lFrame->Gradient->GetMaximum();
1176 float widthFactor = 5.0 + log10f(lFrame->Gradient->GetMaximum());
1177 if (widthFactor < 0.0) widthFactor = 0.0;
1178 VectorScale *= GradMax;
1179 if (Prefs->InvertGradient()) VectorScale *= -1.0;
1180
1181 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, l_specular);
1182 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 30.0f);
1183 Prefs->ChangeColorGradientColor();
1184
1185 CPoint3D NormStart = CPoint3D(0.0f, 0.0f, 1.0f);
1186 for (long iatom=0; iatom<NumAtoms; iatom++) {
1187
1188 if (lAtoms[iatom].GetInvisibility()) continue; //Atom is invisible so skip
1189 CPoint3D NMode;
1190 lFrame->RetrieveAtomGradient(iatom, NMode);
1191 NMode *= VectorScale;
1192 float length = NMode.Magnitude();
1193 if (length > 0.1) {
1194 CPoint3D NModeVector;
1195 Matrix4D rotMat;
1196 NModeVector.x = NMode.x/length;
1197 NModeVector.y = NMode.y/length;
1198 NModeVector.z = NMode.z/length;
1199 float VectorWidth = (0.02+0.03*widthFactor + 0.005*length);
1200 //Set up vectors for the shaft of the normal mode
1201 CPoint3D VStart;
1202 long curAtomType = lAtoms[iatom].Type - 1;
1203
1204 float radius;
1205 if (Prefs->DrawWireFrame()) radius = 0.0;
1206 else radius = AtomScale*Prefs->GetAtomSize(curAtomType);
1207 VStart.x = lAtoms[iatom].Position.x + NModeVector.x*0.95*radius;
1208 VStart.y = lAtoms[iatom].Position.y + NModeVector.y*0.95*radius;
1209 VStart.z = lAtoms[iatom].Position.z + NModeVector.z*0.95*radius;
1210 float HeadRadius = 2 * VectorWidth;
1211 if (2*HeadRadius > length) HeadRadius = length/2.0;
1212 // float HeadRatio = (length-HeadRadius)/length;
1213 GLfloat ShaftLength = length - HeadRadius;
1214
1215 SetRotationMatrix(rotMat, &NormStart, &NModeVector);
1216 rotMat[3][0] = VStart.x;
1217 rotMat[3][1] = VStart.y;
1218 rotMat[3][2] = VStart.z;
1219 glPushMatrix();
1220 glMultMatrixf((const GLfloat *) &rotMat);
1221
1222 gluCylinder(core_obj, VectorWidth, VectorWidth, ShaftLength, (long)(Quality), (long)(0.5*Quality));
1223 glPopMatrix();
1224 rotMat[3][0] = VStart.x + NModeVector.x * ShaftLength;
1225 rotMat[3][1] = VStart.y + NModeVector.y * ShaftLength;
1226 rotMat[3][2] = VStart.z + NModeVector.z * ShaftLength;
1227 glPushMatrix();
1228 glMultMatrixf((const GLfloat *) &rotMat);
1229 gluDisk(core_obj, 0.0, 2*VectorWidth, (long)(Quality), 2);
1230 gluCylinder(core_obj, 2*VectorWidth, 0.0, HeadRadius, (long)(Quality), 3);
1231 glPopMatrix();
1232 }
1233 }
1234 }
1235 gluDeleteQuadric(core_obj); //finally delete the quadric object
1236 }
1237
ChangeColorBondColor(long order) const1238 void WinPrefs::ChangeColorBondColor(long order) const {
1239 float red, green, blue;
1240 red = BondColors[order].red/65536.0;
1241 green = BondColors[order].green/65536.0;
1242 blue = BondColors[order].blue/65536.0;
1243 glColor3f(red, green, blue);
1244 }
1245
ChangeColorAtomColor(long atomtype,float alpha) const1246 void WinPrefs::ChangeColorAtomColor(long atomtype, float alpha) const {
1247 float red, green, blue;
1248 red = AtomColors[atomtype-1].red/65536.0;
1249 green = AtomColors[atomtype-1].green/65536.0;
1250 blue = AtomColors[atomtype-1].blue/65536.0;
1251 glColor4f(red, green, blue, alpha);
1252 }
1253
ChangeColorVectorColor(void) const1254 void WinPrefs::ChangeColorVectorColor(void) const {
1255 float red, green, blue;
1256 red = VectorColor.red/65536.0;
1257 green = VectorColor.green/65536.0;
1258 blue = VectorColor.blue/65536.0;
1259 glColor3f(red, green, blue);
1260 }
1261
ChangeColorGradientColor(void) const1262 void WinPrefs::ChangeColorGradientColor(void) const {
1263 float red, green, blue;
1264 red = GradientColor.red/65536.0;
1265 green = GradientColor.green/65536.0;
1266 blue = GradientColor.blue/65536.0;
1267 glColor3f(red, green, blue);
1268 }
1269
GetAtomColorInverse(long atomtype,float rgb[3]) const1270 void WinPrefs::GetAtomColorInverse(long atomtype, float rgb[3]) const {
1271 rgb[0] = 1.0f - AtomColors[atomtype-1].red / 65536.0f;
1272 rgb[1] = 1.0f - AtomColors[atomtype-1].green / 65536.0f;
1273 rgb[2] = 1.0f - AtomColors[atomtype-1].blue / 65536.0f;
1274 }
1275
DrawHydrogenBond(const Bond & bond,const mpAtom & atom1,const mpAtom & atom2,const WinPrefs & Prefs,GLUquadric *,GLuint sphere_list,bool highlighting_on)1276 void MolDisplayWin::DrawHydrogenBond(const Bond& bond, const mpAtom& atom1,
1277 const mpAtom& atom2,
1278 const WinPrefs& Prefs,
1279 GLUquadric */*quadric*/, GLuint sphere_list,
1280 bool highlighting_on) {
1281 CPoint3D v1, v2, offset;
1282
1283 Prefs.ChangeColorBondColor(kHydrogenBond);
1284 GLdouble BondSize = Prefs.GetQD3DBondWidth() * 0.5;
1285 offset = atom2.Position - atom1.Position;
1286
1287 //Plot as a series of spheres
1288 CPoint3D NormalOffset, NormEnd, NormStart = CPoint3D(0.0f, 0.0f, 1.0f);
1289 Matrix4D rotMat;
1290
1291 float length = offset.Magnitude();
1292 if (length>0.00001) {
1293 NormalOffset.x = offset.x/length;
1294 NormalOffset.y = offset.y/length;
1295 NormalOffset.z = offset.z/length;
1296 } else
1297 NormalOffset.x=NormalOffset.y=NormalOffset.z=0.0;
1298 NormEnd = v2;
1299 Normalize3D(&NormEnd);
1300 SetRotationMatrix(rotMat, &NormStart, &NormalOffset);
1301 rotMat[3][0] = atom1.Position.x;
1302 rotMat[3][1] = atom1.Position.y;
1303 rotMat[3][2] = atom1.Position.z;
1304
1305 glPushMatrix();
1306 glMultMatrixf((const GLfloat *) &rotMat);
1307 float pos=0.75*BondSize;
1308 glTranslatef(0.0, 0.0, pos);
1309 while (pos < length) {
1310 glPushMatrix();
1311 glScalef(BondSize, BondSize, BondSize);
1312 Prefs.ChangeColorBondColor(kHydrogenBond);
1313 glCallList(sphere_list);
1314 if (highlighting_on && !bond.GetSelectState()) {
1315 glColor3f(0.0f, 0.0f, 0.0f);
1316 glEnable(GL_POLYGON_STIPPLE);
1317 glPolygonStipple(stippleMask);
1318
1319 glPushMatrix();
1320 glScalef(1.01f, 1.01f, 1.01f);
1321 glCallList(sphere_list);
1322 glPopMatrix();
1323 glDisable(GL_POLYGON_STIPPLE);
1324 }
1325 glPopMatrix();
1326 glTranslatef(0.0, 0.0, 2.5*BondSize);
1327 pos += 2.5*BondSize;
1328 }
1329 glPopMatrix();
1330 }
1331
AddAxisGL(void)1332 void MolDisplayWin::AddAxisGL(void) {
1333 glEnable(GL_COLOR_MATERIAL);
1334 glColor3fv(fg_color);
1335
1336 glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, l_emissive);
1337 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, l_specular);
1338 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 30.0f);
1339
1340 long Quality = (long)(Prefs->GetQD3DAtomQuality());
1341 float VectorWidth = 0.02;
1342
1343 CPoint3D vector = CPoint3D(1.0f, 0.0f, 0.0f);
1344 CPoint3D NormStart = CPoint3D(0.0f, 0.0f, 1.0f);
1345
1346 // Z-axis
1347 glPushMatrix();
1348 glTranslatef(0.0, 0.0, -MainData->MaxSize);
1349 DrawArrow(2*MainData->MaxSize, VectorWidth, Quality);
1350 glTranslatef(0.0, 0.0, 2.0f * MainData->MaxSize);
1351 DrawSceneString(0.3f, 0.0f, 0.0f, 0.0f, wxT("z"));
1352 glPopMatrix();
1353
1354 // X-axis
1355 Matrix4D rotMat;
1356 SetRotationMatrix(rotMat, &NormStart, &vector);
1357 rotMat[3][0] = -MainData->MaxSize;
1358 glPushMatrix();
1359 glMultMatrixf((const GLfloat *) &rotMat);
1360 DrawArrow(2*MainData->MaxSize, VectorWidth, Quality);
1361 glTranslatef(0.0, 0.0, 2.0f * MainData->MaxSize);
1362 DrawSceneString(0.3f, 0.0f, 0.0f, 0.0f, wxT("x"));
1363 glPopMatrix();
1364
1365 // Y-axis
1366 vector.x = 0.0;
1367 vector.y = 1.0;
1368 SetRotationMatrix(rotMat, &NormStart, &vector);
1369 rotMat[3][1] = -MainData->MaxSize;
1370 glPushMatrix();
1371 glMultMatrixf((const GLfloat *) &rotMat);
1372 DrawArrow(2*MainData->MaxSize, VectorWidth, Quality);
1373 glTranslatef(0.0, 0.0, 2.0f * MainData->MaxSize);
1374 DrawSceneString(0.3f, 0.0f, 0.0f, 0.0f, wxT("y"));
1375 glPopMatrix();
1376 }
1377
AddSymmetryOperators(void)1378 void MolDisplayWin::AddSymmetryOperators(void) {
1379 //Add planes, rotation axis's, etc as appropriate for the selected symmetry point group
1380 //This routine attempts to match the GAMESS assumptions on symmetry
1381 // z is the principle rotation axis (if any)
1382 // x is a perpendicular two-fold axis (if any),
1383 // xz is the sigma-v plane (if any), and
1384 // xy is the sigma-h plane (if any).
1385
1386 if (!MainData->InputOptions) return;
1387 if (!MainData->InputOptions->Data) return;
1388 CPoint3D origin = CPoint3D(0.0f, 0.0f, 0.0f);
1389 CPoint3D p1 = CPoint3D(0.0f, 0.0f, 0.0f), p2 = CPoint3D(0.0f, 0.0f, 0.0f);
1390 switch (MainData->InputOptions->Data->GetPointGroup()) {
1391 case GAMESS_CS: //single sigma-h plane - XY plane
1392 origin.x = origin.y = - MainData->MaxSize;
1393 p1.x = MainData->MaxSize;
1394 p1.y = - p1.x;
1395 p2.x = -MainData->MaxSize;
1396 p2.y = -p2.x;
1397 DrawTranslucentPlane(origin, p1, p2);
1398 break;
1399 case GAMESS_CI:
1400 DrawInversionPoint();
1401 break;
1402 case GAMESS_CNH:
1403 {
1404 int order = MainData->InputOptions->Data->GetPointGroupOrder();
1405 if (order <= 2) order = 2;
1406 if ((order == 2)||(order == 4)||(order == 6)) DrawInversionPoint();
1407 origin.x = origin.y = - MainData->MaxSize;
1408 p1.x = MainData->MaxSize;
1409 p1.y = - p1.x;
1410 p2.x = -MainData->MaxSize;
1411 p2.y = -p2.x;
1412 DrawTranslucentPlane(origin, p1, p2);
1413 origin.x = origin.y = p1.x = p1.y = 0.0;
1414 origin.z = -MainData->MaxSize;
1415 p1.z = MainData->MaxSize;
1416 DrawRotationAxis(origin, p1, order);
1417 }
1418 break;
1419 case GAMESS_CNV: //CN axis plus N sigma-v planes
1420 {
1421 int order = MainData->InputOptions->Data->GetPointGroupOrder();
1422 if (order <= 2) order = 2;
1423 for (int i=0; i<order; i++) {
1424 origin.x = - cos(kPi * i/order) * MainData->MaxSize;
1425 origin.y = - sin(kPi * i/order) * MainData->MaxSize;
1426 origin.z = - MainData->MaxSize;
1427 p1.x = - origin.x;
1428 p1.y = - origin.y;
1429 p1.z = origin.z;
1430 p2.x = origin.x;
1431 p2.y = origin.y;
1432 p2.z = - origin.z;
1433 DrawTranslucentPlane(origin, p1, p2);
1434 }
1435 origin.x = origin.y = p1.x = p1.y = 0.0;
1436 }
1437 // break; //let this roll through to pick up the axis
1438 case GAMESS_CN:
1439 {
1440 int order = MainData->InputOptions->Data->GetPointGroupOrder();
1441 if (order <= 2) order = 2;
1442 origin.z = -MainData->MaxSize;
1443 p1.z = MainData->MaxSize;
1444 DrawRotationAxis(origin, p1, order);
1445 }
1446 break;
1447 case GAMESS_S2N:
1448 break;
1449 case GAMESS_DND:
1450 {
1451 int order = MainData->InputOptions->Data->GetPointGroupOrder();
1452 if (order <= 2) order = 2;
1453 if ((order == 3)||(order == 5)) DrawInversionPoint();
1454 origin.x = origin.y = p1.x = p1.y = 0.0;
1455 origin.z = -MainData->MaxSize;
1456 p1.z = MainData->MaxSize;
1457 DrawRotationAxis(origin, p1, order);
1458 //Set of order C2 axis perpendicular to primary axis
1459 for (int i=0; i<order; i++) {
1460 origin.x = - cos(kPi * i/order) * MainData->MaxSize;
1461 origin.y = - sin(kPi * i/order) * MainData->MaxSize;
1462 origin.z = 0.0;
1463 p1.x = - origin.x;
1464 p1.y = - origin.y;
1465 p1.z = origin.z;
1466 DrawRotationAxis(origin, p1, 2);
1467 }
1468 //There are order count of planes bisecting the C2 axis
1469 for (int i=0; i<order; i++) {
1470 origin.x = - cos(kPi * i/order + kPi/(2*order)) * MainData->MaxSize;
1471 origin.y = - sin(kPi * i/order + kPi/(2*order)) * MainData->MaxSize;
1472 origin.z = -MainData->MaxSize;
1473 p1.x = - origin.x;
1474 p1.y = - origin.y;
1475 p1.z = origin.z;
1476 p2.x = origin.x;
1477 p2.y = origin.y;
1478 p2.z = - origin.z;
1479 DrawTranslucentPlane(origin, p1, p2);
1480 }
1481 }
1482 break;
1483 case GAMESS_DNH:
1484 {
1485 int order = MainData->InputOptions->Data->GetPointGroupOrder();
1486 if (order <= 2) order = 2;
1487 if ((order == 2)||(order == 4)||(order == 6)) DrawInversionPoint();
1488 origin.x = origin.y = - MainData->MaxSize;
1489 p1.x = MainData->MaxSize;
1490 p1.y = - p1.x;
1491 p2.x = -MainData->MaxSize;
1492 p2.y = -p2.x;
1493 DrawTranslucentPlane(origin, p1, p2);
1494 }
1495 // break;
1496 case GAMESS_DN:
1497 {
1498 int order = MainData->InputOptions->Data->GetPointGroupOrder();
1499 if (order <= 2) order = 2;
1500 origin.x = origin.y = p1.x = p1.y = 0.0;
1501 origin.z = -MainData->MaxSize;
1502 p1.z = MainData->MaxSize;
1503 DrawRotationAxis(origin, p1, order);
1504 //Set of order C2 axis perpendicular to primary axis
1505 for (int i=0; i<order; i++) {
1506 origin.x = - cos(kPi * i/order) * MainData->MaxSize;
1507 origin.y = - sin(kPi * i/order) * MainData->MaxSize;
1508 origin.z = 0.0;
1509 p1.x = - origin.x;
1510 p1.y = - origin.y;
1511 p1.z = origin.z;
1512 DrawRotationAxis(origin, p1, 2);
1513 }
1514 }
1515 break;
1516 /* case GAMESS_TD:
1517 break;
1518 case GAMESS_TH:
1519 {
1520 DrawInversionPoint();
1521 }
1522 break;
1523 case GAMESS_T:
1524 break;
1525 case GAMESS_OH:
1526 {
1527 DrawInversionPoint();
1528 }
1529 break;
1530 case GAMESS_O:
1531 break;
1532 */
1533 default: //The rest are not implemented.
1534 break;
1535 }
1536 }
1537
Draw3DGL(MoleculeData * MainData,WinPrefs *,myGLTriangle *,unsigned int)1538 long Surf1DBase::Draw3DGL(MoleculeData *MainData, WinPrefs */*Prefs*/, myGLTriangle *, unsigned int /*shader_program*/) {
1539
1540 if (Visible) {
1541 // Update the grid if needed, then display
1542 if (!Grid) Update(MainData);
1543 }
1544
1545 glLoadName(MMP_1DSURFACE);
1546 glPushName(GetSurfaceID()+1);
1547
1548 glPushName(1);
1549
1550 glPushMatrix();
1551 glTranslatef(Start.x, Start.y, Start.z);
1552 glScalef(0.05f, 0.05f, 0.05f);
1553 glCallList(MainData->MolWin->GetSphereList());
1554 glPopMatrix();
1555
1556 glLoadName(2);
1557 glPushMatrix();
1558 glTranslatef(End.x, End.y, End.z);
1559 glScalef(0.05f, 0.05f, 0.05f);
1560 glCallList(MainData->MolWin->GetSphereList());
1561 glPopMatrix();
1562
1563 CPoint3D vec = End - Start;
1564 CPoint3D perp;
1565 OrthoVector(vec, perp);
1566 vec *= 1.0f / (NumGridPoints - 1);
1567 CPoint3D curr_pt = Start;
1568 CPoint3D next_pt;
1569
1570 // Find a vector perpendicular to both view direction and the vector from
1571 // Start to End projected onto the screen plane. The graph will be pushed
1572 // away along this perpendicular vector so that it can always be seen, no
1573 // matter the rotation.
1574 const Matrix4D& rot = MainData->GetRotationMatrix();
1575 CPoint3D normal(rot[0][2], rot[1][2], rot[2][2]);
1576 CPoint3D projected_start = Start;
1577 CPoint3D projected_end = End;
1578 ProjectToPlane(normal, Start, projected_start);
1579 ProjectToPlane(normal, Start, projected_end);
1580 CPoint3D projected_vec = projected_end - projected_start;
1581 CrossProduct3D(&projected_vec, &normal, &perp);
1582 Normalize3D(&perp);
1583
1584 glLoadName(0);
1585 float range_factor;
1586 if (GridMax > MaxContourValue) {
1587 range_factor = 1.0f / (MaxContourValue - GridMin);
1588 } else {
1589 range_factor = 1.0f / (GridMax - GridMin);
1590 }
1591
1592 float red;
1593 float curr_val, next_val;
1594 curr_val = Grid[0] > MaxContourValue ? MaxContourValue : Grid[0];
1595
1596 #if 1
1597 GLUquadricObj *quadric = gluNewQuadric();
1598 CreateCylinderFromLine(quadric, Start, End, 0.02f, 5, 5, false);
1599
1600 for (int i = 0; i < NumGridPoints - 1; ++i) {
1601 next_val = Grid[i + 1] > MaxContourValue ? MaxContourValue : Grid[i + 1];
1602 next_pt = curr_pt + vec;
1603 red = (curr_val - GridMin) * range_factor;
1604 glColor3f(red, 0.0f, 1.0f - red);
1605 CreateCylinderFromLine(quadric, curr_pt + perp * (curr_val * Scale),
1606 next_pt + perp * (next_val * Scale),
1607 0.05f, 5, 5, true);
1608 curr_val = next_val;
1609 curr_pt = next_pt;
1610 }
1611
1612 glColor3f(1.0f, 1.0f, 1.0f);
1613 if (fabs(Grid[0]) > 1e-6f) {
1614 glPushMatrix();
1615 glTranslatef(next_pt.x, next_pt.y, next_pt.z);
1616 glScalef(0.05f, 0.05f, 0.05f);
1617 glCallList(MainData->MolWin->GetSphereList());
1618 glPopMatrix();
1619
1620 curr_val = Scale * (Grid[0] > MaxContourValue ? MaxContourValue : Grid[0]);
1621 CreateCylinderFromLine(quadric, Start, Start + perp * curr_val,
1622 0.02f, 5, 5, true);
1623 }
1624
1625 if (fabs(Grid[NumGridPoints - 1]) > 1e-6f) {
1626 curr_val = Scale * (Grid[NumGridPoints - 1] > MaxContourValue ? MaxContourValue : Grid[NumGridPoints - 1]);
1627 CreateCylinderFromLine(quadric, End, End + perp * curr_val,
1628 0.02f, 5, 5, true);
1629 }
1630
1631 gluDeleteQuadric(quadric);
1632 #else
1633 if (GLEW_VERSION_2_0 && Prefs->GetShaderMode()) {
1634 glUseProgram(0);
1635 }
1636
1637 glDisable(GL_LIGHTING);
1638 glBegin(GL_QUADS);
1639 for (int i = 0; i < NumGridPoints - 1; ++i) {
1640 next_val = Grid[i + 1] > MaxContourValue ? MaxContourValue : Grid[i + 1];
1641 next_pt = curr_pt + vec;
1642 red = (curr_val - GridMin) * range_factor;
1643 glColor3f(red, 0.0f, 1.0f - red);
1644 glVertex3f(curr_pt.x - perp.x * 0.045f,
1645 curr_pt.y - perp.y * 0.045f,
1646 curr_pt.z - perp.z * 0.045f);
1647 glVertex3f(curr_pt.x + perp.x * 0.045f,
1648 curr_pt.y + perp.y * 0.045f,
1649 curr_pt.z + perp.z * 0.045f);
1650 red = (Grid[i + 1] - GridMin) * range_factor;
1651 glColor3f(red, 0.0f, 1.0f - red);
1652 glVertex3f(next_pt.x + perp.x * 0.045f,
1653 next_pt.y + perp.y * 0.045f,
1654 next_pt.z + perp.z * 0.045f);
1655 glVertex3f(next_pt.x - perp.x * 0.045f,
1656 next_pt.y - perp.y * 0.045f,
1657 next_pt.z - perp.z * 0.045f);
1658 curr_pt = next_pt;
1659 curr_val = next_val;
1660 }
1661 glEnd();
1662 glEnable(GL_LIGHTING);
1663
1664 if (GLEW_VERSION_2_0 && Prefs->GetShaderMode()) {
1665 glUseProgram(shader_program);
1666 }
1667 #endif
1668
1669 glPopName();
1670 glPopName();
1671 return 0;
1672 }
1673
Draw3DGL(MoleculeData * MainData,WinPrefs * Prefs,myGLTriangle *,unsigned int)1674 long Surf2DBase::Draw3DGL(MoleculeData * MainData, WinPrefs * Prefs, myGLTriangle *, unsigned int /*shader_program*/) {
1675 if (Visible) {
1676 //Update the grid if needed, then contour and display
1677 if (!Grid) Update(MainData);
1678 if (Grid) Contour2DGrid3DGL(MainData, Prefs);
1679 if (ShowPlottingPlane()) DrawTranslucentPlane(Origin, Origin+XInc*(NumGridPoints-1), Origin+YInc*(NumGridPoints-1));
1680 }
1681 return 0;
1682 }
1683
Contour2DGrid3DGL(MoleculeData *,WinPrefs * Prefs)1684 void Surf2DBase::Contour2DGrid3DGL(MoleculeData * , WinPrefs * Prefs)
1685 {
1686 //Scan the Grid producing the contours
1687 float TestPoint1, TestPoint2, TestPoint3, TestPoint4, XGridValue, YGridValue, ZGridValue;
1688 CPoint3D Contour[4];
1689 Boolean HasPoint[4];
1690
1691 long NumPoints = NumGridPoints;
1692 float * lGrid = Grid;
1693
1694 CPoint3D XGridMin, XGridInc, YGridInc;
1695 XGridMin = Origin;
1696 // XGridMin *= kBohr2AngConversion;
1697 XGridInc = XInc;
1698 // XGridInc *= kBohr2AngConversion;
1699 YGridInc = YInc;
1700 // YGridInc *= kBohr2AngConversion;
1701
1702 float ContourValueInc = MaxContourValue/(NumContours+1);
1703 long NumPosContours = (long)(fabs(GridMax)/ContourValueInc) + 1;
1704 long NumNegContours = (long)(fabs(GridMin)/ContourValueInc) + 1;
1705
1706 CPoint3D lineStart, lineEnd;
1707 float ContourValue = 0.0;
1708 float lineWidth = Prefs->GetQD3DLineWidth();
1709 GLUquadricObj * qobj = NULL;
1710 Boolean UseLines = true;
1711 if (lineWidth > 0.0001) {
1712 UseLines = false;
1713 qobj = gluNewQuadric();
1714 if (!qobj) throw std::bad_alloc();
1715 }
1716 long n;
1717 Boolean Dash = GetDashLine();
1718 // GLenum error = glGetError();
1719 if (UseLines) {
1720 glDisable(GL_LIGHTING);
1721 glLineWidth(1);
1722 glBegin(GL_LINES);
1723 // error = glGetError(); //maybe clear off odd errors...
1724 }
1725 glColor3f(0.65, 0.65, 0.65); //Setup for the gray zero contour color
1726
1727 //Go up to NumContours+1 to allow for the zero contour
1728 for (long iContour=0; iContour<=NumContours; iContour++) {
1729 for (int pass=0; pass<2; pass++) {
1730 if (iContour==0) { //0 value contour
1731 pass++; //only need one 0 contour!
1732 //Plot zero value contour only if requested
1733 if (!(SurfOptions&1)) continue;
1734 } else {
1735 ContourValue *= -1.0;
1736 if (pass==0) {
1737 if (!ContourBothPosNeg()) continue;
1738 if (iContour > NumNegContours) continue;
1739 glColor3f((float)NegColor.red/65536, (float)NegColor.green/65536, (float)NegColor.blue/65536);
1740 } else {
1741 if (iContour > NumPosContours) continue;
1742 glColor3f((float)PosColor.red/65536, (float)PosColor.green/65536, (float)PosColor.blue/65536);
1743 }
1744 }
1745 n=NumGridPoints;
1746 for (long i=1; i<NumPoints; i++) {
1747 XGridValue = XGridMin.x + i*XGridInc.x;
1748 YGridValue = XGridMin.y + i*XGridInc.y;
1749 ZGridValue = XGridMin.z + i*XGridInc.z;
1750 for (long j=1; j<NumPoints; j++) {
1751 XGridValue += YGridInc.x;
1752 YGridValue += YGridInc.y;
1753 ZGridValue += YGridInc.z;
1754 n++;
1755 for (int i=0; i<4; i++) HasPoint[i]=false;
1756
1757 TestPoint1 = lGrid[n]-ContourValue;
1758 TestPoint2 = lGrid[n-1]-ContourValue;
1759 TestPoint3 = lGrid[n-NumPoints]-ContourValue;
1760 TestPoint4 = lGrid[n-1-NumPoints]-ContourValue;
1761
1762 if ((TestPoint1*TestPoint2)<0.0) {
1763 HasPoint[0]=true;
1764 Contour[0].x = XGridValue - YGridInc.x*(TestPoint1/(TestPoint1-TestPoint2));
1765 Contour[0].y = YGridValue - YGridInc.y*(TestPoint1/(TestPoint1-TestPoint2));
1766 Contour[0].z = ZGridValue - YGridInc.z*(TestPoint1/(TestPoint1-TestPoint2));
1767 }
1768 if ((TestPoint1*TestPoint3)<0.0) {
1769 HasPoint[1]=true;
1770 Contour[1].x = XGridValue - XGridInc.x*(TestPoint1/(TestPoint1-TestPoint3));
1771 Contour[1].y = YGridValue - XGridInc.y*(TestPoint1/(TestPoint1-TestPoint3));
1772 Contour[1].z = ZGridValue - XGridInc.z*(TestPoint1/(TestPoint1-TestPoint3));
1773 }
1774 if ((TestPoint2*TestPoint4)<0.0) {
1775 HasPoint[2]=true;
1776 Contour[2].x = XGridValue-YGridInc.x - XGridInc.x*(TestPoint2/(TestPoint2-TestPoint4));
1777 Contour[2].y = YGridValue-YGridInc.y - XGridInc.y*(TestPoint2/(TestPoint2-TestPoint4));
1778 Contour[2].z = ZGridValue-YGridInc.z - XGridInc.z*(TestPoint2/(TestPoint2-TestPoint4));
1779 }
1780 if ((TestPoint3*TestPoint4)<0.0) {
1781 HasPoint[3]=true;
1782 Contour[3].x = XGridValue-XGridInc.x - YGridInc.x*(TestPoint3/(TestPoint3-TestPoint4));
1783 Contour[3].y = YGridValue-XGridInc.y - YGridInc.y*(TestPoint3/(TestPoint3-TestPoint4));
1784 Contour[3].z = ZGridValue-XGridInc.z - YGridInc.z*(TestPoint3/(TestPoint3-TestPoint4));
1785 }
1786
1787 if (HasPoint[0]) {
1788 lineStart = Contour[0];
1789 if (HasPoint[1] || HasPoint[2]) {
1790 if (HasPoint[1]) {
1791 lineEnd = Contour[1];
1792 if ((pass==0)&&Dash) {
1793 lineEnd.x = lineStart.x + (lineEnd.x-lineStart.x)/2.0;
1794 lineEnd.y = lineStart.y + (lineEnd.y-lineStart.y)/2.0;
1795 lineEnd.z = lineStart.z + (lineEnd.z-lineStart.z)/2.0;
1796 }
1797 if (UseLines) {
1798 glVertex3d(lineStart.x, lineStart.y, lineStart.z);
1799 glVertex3d(lineEnd.x, lineEnd.y, lineEnd.z);
1800 } else {
1801 CreateCylinderFromLine(qobj, lineStart, lineEnd, lineWidth);
1802 }
1803 }
1804 if (HasPoint[2]) {
1805 lineEnd = Contour[2];
1806 if ((pass==0)&&Dash) {
1807 lineEnd.x = lineStart.x + (lineEnd.x-lineStart.x)/2.0;
1808 lineEnd.y = lineStart.y + (lineEnd.y-lineStart.y)/2.0;
1809 lineEnd.z = lineStart.z + (lineEnd.z-lineStart.z)/2.0;
1810 }
1811 if (UseLines) {
1812 glVertex3d(lineStart.x, lineStart.y, lineStart.z);
1813 glVertex3d(lineEnd.x, lineEnd.y, lineEnd.z);
1814 } else {
1815 CreateCylinderFromLine(qobj, lineStart, lineEnd, lineWidth);
1816 }
1817 }
1818 } else if (HasPoint[3]) {
1819 lineEnd = Contour[3];
1820 if ((pass==0)&&Dash) {
1821 lineEnd.x = lineStart.x + (lineEnd.x-lineStart.x)/2.0;
1822 lineEnd.y = lineStart.y + (lineEnd.y-lineStart.y)/2.0;
1823 lineEnd.z = lineStart.z + (lineEnd.z-lineStart.z)/2.0;
1824 }
1825 if (UseLines) {
1826 glVertex3d(lineStart.x, lineStart.y, lineStart.z);
1827 glVertex3d(lineEnd.x, lineEnd.y, lineEnd.z);
1828 } else {
1829 CreateCylinderFromLine(qobj, lineStart, lineEnd, lineWidth);
1830 }
1831 }
1832 }
1833 if (HasPoint[1]) {
1834 if (HasPoint[3]) {
1835 lineStart = Contour[1];
1836 lineEnd = Contour[3];
1837 if ((pass==0)&&Dash) {
1838 lineEnd.x = lineStart.x + (lineEnd.x-lineStart.x)/2.0;
1839 lineEnd.y = lineStart.y + (lineEnd.y-lineStart.y)/2.0;
1840 lineEnd.z = lineStart.z + (lineEnd.z-lineStart.z)/2.0;
1841 }
1842 if (UseLines) {
1843 glVertex3d(lineStart.x, lineStart.y, lineStart.z);
1844 glVertex3d(lineEnd.x, lineEnd.y, lineEnd.z);
1845 } else {
1846 CreateCylinderFromLine(qobj, lineStart, lineEnd, lineWidth);
1847 }
1848 } else if (HasPoint[2] && !HasPoint[0]) {
1849 lineStart = Contour[1];
1850 lineEnd = Contour[2];
1851 if ((pass==0)&&Dash) {
1852 lineEnd.x = lineStart.x + (lineEnd.x-lineStart.x)/2.0;
1853 lineEnd.y = lineStart.y + (lineEnd.y-lineStart.y)/2.0;
1854 lineEnd.z = lineStart.z + (lineEnd.z-lineStart.z)/2.0;
1855 }
1856 if (UseLines) {
1857 glVertex3d(lineStart.x, lineStart.y, lineStart.z);
1858 glVertex3d(lineEnd.x, lineEnd.y, lineEnd.z);
1859 } else {
1860 CreateCylinderFromLine(qobj, lineStart, lineEnd, lineWidth);
1861 }
1862 }
1863 }
1864 if (HasPoint[2]&&HasPoint[3]) {
1865 lineStart = Contour[2];
1866 lineEnd = Contour[3];
1867 if ((pass==0)&&Dash) {
1868 lineEnd.x = lineStart.x + (lineEnd.x-lineStart.x)/2.0;
1869 lineEnd.y = lineStart.y + (lineEnd.y-lineStart.y)/2.0;
1870 lineEnd.z = lineStart.z + (lineEnd.z-lineStart.z)/2.0;
1871 }
1872 if (UseLines) {
1873 glVertex3d(lineStart.x, lineStart.y, lineStart.z);
1874 glVertex3d(lineEnd.x, lineEnd.y, lineEnd.z);
1875 } else {
1876 CreateCylinderFromLine(qobj, lineStart, lineEnd, lineWidth);
1877 }
1878 }
1879 }
1880 n++;
1881 }
1882 }
1883 ContourValue += ContourValueInc;
1884 }
1885 if (UseLines) {
1886 glEnd();
1887 // error = glGetError(); //This is here to clear off odd errors.
1888 glEnable(GL_LIGHTING);
1889 } else {
1890 if (qobj) gluDeleteQuadric(qobj); //finally delete the quadric object
1891 }
1892 }
1893
getTriangleCount(void) const1894 long General3DSurface::getTriangleCount(void) const {
1895 long result = NumPosContourTriangles;
1896 if (Mode & 4) result += NumNegContourTriangles;
1897 return result;
1898 }
1899
Draw3DGL(MoleculeData * MainData,WinPrefs * Prefs,myGLTriangle * transpTri,unsigned int)1900 long General3DSurface::Draw3DGL(MoleculeData * MainData, WinPrefs * Prefs, myGLTriangle * transpTri, unsigned int /*shader_program*/) {
1901 long result=0;
1902 if (Visible) {
1903 if (ContourHndl && VertexList) {
1904 if (SolidSurface()) {
1905 if (UseSurfaceNormals()&&SurfaceNormals) {
1906 result = CreateSolidSurface(ContourHndl, SurfaceNormals, VertexList,
1907 NumPosContourTriangles, &PosColor, List, &NegColor, MaxMEPValue, MainData, transpTri);
1908 if ((Mode & 4)&&(NumNegContourTriangles > 0))
1909 result += CreateSolidSurface(ContourHndl,
1910 SurfaceNormals, &(VertexList[3*NumPosContourTriangles]),
1911 NumNegContourTriangles, &NegColor, List, &PosColor, MaxMEPValue, MainData, &(transpTri[result]));
1912 } else {
1913 result = CreateSolidSurface(ContourHndl, NULL, VertexList,
1914 NumPosContourTriangles, &PosColor, List, &NegColor, MaxMEPValue, MainData, transpTri);
1915 if ((Mode & 4)&&(NumNegContourTriangles > 0))
1916 result += CreateSolidSurface(ContourHndl,
1917 NULL, &(VertexList[3*NumPosContourTriangles]),
1918 NumNegContourTriangles, &NegColor, List, &PosColor, MaxMEPValue, MainData, &(transpTri[result]));
1919 }
1920 } else if (WireFrameSurface()) {
1921 CreateWireSurface(ContourHndl, NULL, VertexList,
1922 NumPosContourTriangles, &PosColor, List, &NegColor, MaxMEPValue, MainData, Prefs);
1923 if (ContourBothPosNeg()&&(NumNegContourTriangles > 0))
1924 CreateWireSurface(ContourHndl, NULL,
1925 &(VertexList[3*NumPosContourTriangles]),
1926 NumNegContourTriangles, &NegColor, List, &PosColor, MaxMEPValue, MainData, Prefs);
1927 }
1928 }
1929 }
1930 return result;
1931 }
Draw3DGL(MoleculeData * MainData,WinPrefs * Prefs,myGLTriangle * transpTri,unsigned int)1932 long TEDensity3DSurface::Draw3DGL(MoleculeData * MainData, WinPrefs * Prefs, myGLTriangle * transpTri, unsigned int /*shader_program*/) {
1933 long result = 0;
1934 if (Visible) {
1935 if (ContourHndl && VertexList) {
1936 if (SolidSurface()) {
1937 if ((UseSurfaceNormals())&&SurfaceNormals) {
1938 result = CreateSolidSurface(ContourHndl, SurfaceNormals, VertexList,
1939 NumPosContourTriangles,
1940 &PosColor, List, &NegColor, MaxMEPValue, MainData, transpTri);
1941 } else {
1942 result = CreateSolidSurface(ContourHndl, NULL, VertexList,
1943 NumPosContourTriangles,
1944 &PosColor, List, &NegColor, MaxMEPValue, MainData, transpTri);
1945 }
1946 } else if (WireFrameSurface()) {
1947 CreateWireSurface(ContourHndl, NULL, VertexList,
1948 NumPosContourTriangles,
1949 &PosColor, List, &NegColor, MaxMEPValue, MainData, Prefs);
1950 }
1951 }
1952 }
1953 return result;
1954 }
Draw3DGL(MoleculeData * MainData,WinPrefs * Prefs,myGLTriangle * transpTri,unsigned int)1955 long Orb3DSurface::Draw3DGL(MoleculeData * MainData, WinPrefs * Prefs, myGLTriangle * transpTri, unsigned int /*shader_program*/) {
1956 long result=0;
1957 if (Visible && (PlotOrb>=0)) {
1958 if (ContourHndl && VertexList) {
1959 if (SolidSurface()) {
1960 if ((UseSurfaceNormals())&&SurfaceNormals) {
1961 result = CreateSolidSurface(ContourHndl,
1962 SurfaceNormals, VertexList, NumPosContourTriangles,
1963 &PosColor, NULL, NULL, 1.0, MainData, transpTri);
1964 result += CreateSolidSurface(ContourHndl,
1965 SurfaceNormals,
1966 &(VertexList[3*NumPosContourTriangles]),
1967 NumNegContourTriangles, &NegColor,
1968 NULL, NULL, 1.0, MainData, &(transpTri[result]));
1969 } else {
1970 result = CreateSolidSurface(ContourHndl, NULL, VertexList,
1971 NumPosContourTriangles,
1972 &PosColor, NULL, NULL, 1.0, MainData, transpTri);
1973 result += CreateSolidSurface(ContourHndl,
1974 NULL, &(VertexList[3*NumPosContourTriangles]),
1975 NumNegContourTriangles, &NegColor, NULL, NULL, 1.0, MainData, &(transpTri[result]));
1976 }
1977 } else if (WireFrameSurface()) {
1978 CreateWireSurface(ContourHndl, NULL, VertexList, NumPosContourTriangles,
1979 &PosColor, NULL, NULL, 1.0, MainData, Prefs);
1980 CreateWireSurface(ContourHndl, NULL,
1981 &(VertexList[3*NumPosContourTriangles]),
1982 NumNegContourTriangles, &NegColor, NULL, NULL, 1.0, MainData, Prefs);
1983 }
1984 }
1985 }
1986 return result;
1987 }
Draw3DGL(MoleculeData * MainData,WinPrefs * Prefs,myGLTriangle * transpTri,unsigned int)1988 long MEP3DSurface::Draw3DGL(MoleculeData * MainData, WinPrefs * Prefs, myGLTriangle * transpTri, unsigned int /*shader_program*/) {
1989 long result=0;
1990 if (Visible) {
1991 if (ContourHndl && VertexList) {
1992 if (SolidSurface()) {
1993 if ((UseSurfaceNormals())&&SurfaceNormals) {
1994 result = CreateSolidSurface(ContourHndl, SurfaceNormals,
1995 VertexList, NumPosContourTriangles,
1996 &PosColor, NULL, NULL, 1.0, MainData, transpTri);
1997 result += CreateSolidSurface(ContourHndl, SurfaceNormals,
1998 &(VertexList[3*NumPosContourTriangles]),
1999 NumNegContourTriangles, &NegColor, NULL, NULL, 1.0, MainData, &(transpTri[result]));
2000 } else {
2001 result = CreateSolidSurface(ContourHndl, NULL, VertexList, NumPosContourTriangles,
2002 &PosColor, NULL, NULL, 1.0, MainData, transpTri);
2003 result += CreateSolidSurface(ContourHndl,
2004 NULL, &(VertexList[3*NumPosContourTriangles]),
2005 NumNegContourTriangles, &NegColor, NULL, NULL, 1.0, MainData, &(transpTri[result]));
2006 }
2007 } else if (WireFrameSurface()) {
2008 CreateWireSurface(ContourHndl, NULL, VertexList, NumPosContourTriangles,
2009 &PosColor, NULL, NULL, 1.0, MainData, Prefs);
2010 CreateWireSurface(ContourHndl, NULL,
2011 &(VertexList[3*NumPosContourTriangles]),
2012 NumNegContourTriangles, &NegColor, NULL, NULL, 1.0, MainData, Prefs);
2013 }
2014 }
2015 }
2016 return result;
2017 }
2018
CreateWireSurface(CPoint3D * Vertices,CPoint3D * Normals,long * VertexList,long NumTriangles,RGBColor * SurfaceColor,float * SurfaceValue,RGBColor * NColor,float MaxSurfaceValue,MoleculeData * MainData,const WinPrefs * Prefs)2019 void Surf3DBase::CreateWireSurface(CPoint3D * Vertices, CPoint3D * Normals, long * VertexList,
2020 long NumTriangles, RGBColor * SurfaceColor, float * SurfaceValue,
2021 RGBColor * NColor, float MaxSurfaceValue, MoleculeData * MainData, const WinPrefs * Prefs)
2022 {
2023 glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
2024 glDisable(GL_LIGHTING);
2025 glLineWidth(Prefs->GetLineWidth());
2026 CreateSolidSurface(Vertices, Normals, VertexList,
2027 NumTriangles, SurfaceColor, SurfaceValue, NColor, MaxSurfaceValue, MainData, NULL);
2028 glEnable(GL_LIGHTING);
2029 glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
2030
2031 /* glDisable(GL_LIGHTING);
2032 glLineWidth(1);
2033 //This code works only slightly better than the above code on the Intel Macs with ATI hardware.
2034 //It seems to miss some lines?
2035 CreateWireFrameSurfaceWithLines(Vertices, VertexList,
2036 NumTriangles, SurfaceColor, SurfaceValue, NColor, MaxSurfaceValue, MainData);
2037 // glEnable(GL_LIGHTING);
2038 */
2039 }
2040
SetSurfaceColor(const float & surfaceValue,const RGBColor * pColor,const RGBColor * nColor,float & red,float & green,float & blue) const2041 void Surf3DBase::SetSurfaceColor(const float & surfaceValue, const RGBColor * pColor, const RGBColor * nColor,
2042 float & red, float & green, float & blue) const {
2043 float localVal = surfaceValue;
2044 if (UseRGBColoration()) {
2045 if (InvertRGBColoration()) localVal *= -1.0;
2046 if (localVal < 0.0) {
2047 if (localVal < -1.0) localVal = -1.0;
2048 red = 0.0;
2049 green = localVal + 1.0;
2050 blue = -localVal;
2051 } else {
2052 if (localVal > 1.0) localVal = 1.0;
2053 red = localVal;
2054 green = 1 - localVal;
2055 blue = 0.0;
2056 }
2057 } else { //color based on +/- color intensity
2058 if (localVal>=0) {
2059 if (localVal > 1.0) localVal = 1.0;
2060 if (pColor) {
2061 red = ((float) pColor->red/65536)*localVal;
2062 green = ((float) pColor->green/65536)*localVal;
2063 blue = ((float) pColor->blue/65536)*localVal;
2064 } else {
2065 red = ((float) PosColor.red/65536)*localVal;
2066 green = ((float) PosColor.green/65536)*localVal;
2067 blue = ((float) PosColor.blue/65536)*localVal;
2068 }
2069 } else {
2070 localVal *= -1.0;
2071 if (localVal > 1.0) localVal = 1.0;
2072 if (nColor) {
2073 red = ((float) nColor->red/65536)*localVal;
2074 green = ((float) nColor->green/65536)*localVal;
2075 blue = ((float) nColor->blue/65536)*localVal;
2076 } else { //Hopefully the following is never used
2077 red = ((float) NegColor.red/65536)*localVal;
2078 green = ((float) NegColor.green/65536)*localVal;
2079 blue = ((float) NegColor.blue/65536)*localVal;
2080 }
2081 }
2082 }
2083 }
2084
CreateSolidSurface(CPoint3D * Vertices,CPoint3D * Normals,long * vList,long NumTriangles,RGBColor * SurfaceColor,float * SurfaceValue,RGBColor * NColor,float MaxSurfaceValue,MoleculeData *,myGLTriangle * transpTri)2085 long Surf3DBase::CreateSolidSurface(CPoint3D * Vertices, CPoint3D * Normals, long * vList,
2086 long NumTriangles, RGBColor * SurfaceColor, float * SurfaceValue,
2087 RGBColor * NColor, float MaxSurfaceValue, MoleculeData * , myGLTriangle * transpTri)
2088 {
2089 long v1, v2, v3, result=0;
2090 GLfloat alpha=1.0, red, green, blue, xnorm, ynorm, znorm;
2091
2092 red = (float) SurfaceColor->red/65536.0;
2093 green = (float) SurfaceColor->green/65536.0;
2094 blue = (float) SurfaceColor->blue/65536.0;
2095 red = MIN(red, 1.0);
2096 blue = MIN(blue, 1.0);
2097 green = MIN(green, 1.0);
2098 red = MAX(red, 0.0);
2099 blue = MAX(blue, 0.0);
2100 green = MAX(green, 0.0);
2101 long * VertexList = vList;
2102 if (isTransparent()) {
2103 alpha = (100 - Transparency) / 100.0f;
2104 if (!transpTri) return 0; //transparncy requires a different draw method
2105 result = NumTriangles;
2106 }
2107
2108 if (!SurfaceValue) { //If we are not using surface coloring setup the color once for all the triangles
2109 glColor4f(red, green, blue, alpha);
2110 }
2111
2112
2113 if (!transpTri)
2114 glBegin(GL_TRIANGLES);
2115 for (long itri=0; itri<NumTriangles; itri++) {
2116 v1 = VertexList[3*itri];
2117 v2 = VertexList[3*itri+1];
2118 v3 = VertexList[3*itri+2];
2119
2120 if (Normals) {
2121 xnorm = Normals[v1].x;
2122 ynorm = Normals[v1].y;
2123 znorm = Normals[v1].z;
2124 } else { //compute a simple triangle normal for all three vertices
2125 float qx = Vertices[v2].x - Vertices[v1].x;
2126 float qy = Vertices[v2].y - Vertices[v1].y;
2127 float qz = Vertices[v2].z - Vertices[v1].z;
2128 float px = Vertices[v3].x - Vertices[v1].x;
2129 float py = Vertices[v3].y - Vertices[v1].y;
2130 float pz = Vertices[v3].z - Vertices[v1].z;
2131 xnorm = -(py*qz - pz*qy);
2132 ynorm = -(pz*qx - px*qz);
2133 znorm = -(px*qy - py*qx);
2134
2135 float len = 1.0/sqrt( xnorm*xnorm + ynorm*ynorm + znorm*znorm );
2136 xnorm *= len;
2137 ynorm *= len;
2138 znorm *= len;
2139 }
2140 if (SurfaceValue) {
2141 float temp = SurfaceValue[v1];
2142 temp /= MaxSurfaceValue;
2143 SetSurfaceColor(temp, SurfaceColor, NColor, red, green, blue);
2144 if (!transpTri)
2145 glColor4f(red, green, blue, alpha);
2146 }
2147
2148 if (!transpTri) {
2149 glNormal3f(xnorm, ynorm, znorm);
2150 glVertex3d(Vertices[v1].x, Vertices[v1].y, Vertices[v1].z);
2151 } else {
2152 transpTri[itri].v1 = Vertices[v1];
2153 transpTri[itri].n1.x = xnorm;
2154 transpTri[itri].n1.y = ynorm;
2155 transpTri[itri].n1.z = znorm;
2156 transpTri[itri].r1 = red;
2157 transpTri[itri].g1 = green;
2158 transpTri[itri].b1 = blue;
2159 transpTri[itri].a1 = alpha;
2160 }
2161
2162 if (Normals) {
2163 xnorm = Normals[v2].x;
2164 ynorm = Normals[v2].y;
2165 znorm = Normals[v2].z;
2166 }
2167 if (SurfaceValue) {
2168 float temp = SurfaceValue[v2];
2169 temp /= MaxSurfaceValue;
2170 SetSurfaceColor(temp, SurfaceColor, NColor, red, green, blue);
2171 if (!transpTri)
2172 glColor4f(red, green, blue, alpha);
2173 }
2174 if (!transpTri) {
2175 glNormal3f(xnorm, ynorm, znorm);
2176 glVertex3d(Vertices[v2].x, Vertices[v2].y, Vertices[v2].z);
2177 } else {
2178 transpTri[itri].v2 = Vertices[v2];
2179 transpTri[itri].n2.x = xnorm;
2180 transpTri[itri].n2.y = ynorm;
2181 transpTri[itri].n2.z = znorm;
2182 transpTri[itri].r2 = red;
2183 transpTri[itri].g2 = green;
2184 transpTri[itri].b2 = blue;
2185 transpTri[itri].a2 = alpha;
2186 }
2187
2188 if (Normals) {
2189 xnorm = Normals[v3].x;
2190 ynorm = Normals[v3].y;
2191 znorm = Normals[v3].z;
2192 }
2193 if (SurfaceValue) {
2194 float temp = SurfaceValue[v3];
2195 temp /= MaxSurfaceValue;
2196 SetSurfaceColor(temp, SurfaceColor, NColor, red, green, blue);
2197 if (!transpTri)
2198 glColor4f(red, green, blue, alpha);
2199 }
2200 if (!transpTri) {
2201 glNormal3f(xnorm, ynorm, znorm);
2202 glVertex3d(Vertices[v3].x, Vertices[v3].y, Vertices[v3].z);
2203 } else {
2204 transpTri[itri].v3 = Vertices[v3];
2205 transpTri[itri].n3.x = xnorm;
2206 transpTri[itri].n3.y = ynorm;
2207 transpTri[itri].n3.z = znorm;
2208 transpTri[itri].r3 = red;
2209 transpTri[itri].g3 = green;
2210 transpTri[itri].b3 = blue;
2211 transpTri[itri].a3 = alpha;
2212 }
2213 }
2214 if (!transpTri)
2215 glEnd(); //End of triangle creation
2216 return result;
2217 }
CreateWireFrameSurfaceWithLines(CPoint3D * Vertices,long * vList,long NumTriangles,RGBColor * SurfaceColor,float * SurfaceValue,RGBColor * NColor,float MaxSurfaceValue,MoleculeData *)2218 long Surf3DBase::CreateWireFrameSurfaceWithLines(CPoint3D * Vertices, long * vList,
2219 long NumTriangles, RGBColor * SurfaceColor, float * SurfaceValue,
2220 RGBColor * NColor, float MaxSurfaceValue, MoleculeData * )
2221 {
2222 long v1, v2, v3, result=0;
2223 GLfloat alpha=1.0, red, green, blue;
2224
2225 red = (float) SurfaceColor->red/65536.0;
2226 green = (float) SurfaceColor->green/65536.0;
2227 blue = (float) SurfaceColor->blue/65536.0;
2228 red = MIN(red, 1.0);
2229 blue = MIN(blue, 1.0);
2230 green = MIN(green, 1.0);
2231 red = MAX(red, 0.0);
2232 blue = MAX(blue, 0.0);
2233 green = MAX(green, 0.0);
2234 long * VertexList = vList;
2235
2236 if (!SurfaceValue) { //If we are not using surface coloring setup the color once for all the triangles
2237 glColor4f(red, green, blue, alpha);
2238 }
2239
2240 glLineWidth(1);
2241 glBegin(GL_LINES);
2242 for (long itri=0; itri<NumTriangles; itri++) {
2243 v1 = VertexList[3*itri];
2244 v2 = VertexList[3*itri+1];
2245 v3 = VertexList[3*itri+2];
2246
2247 if (SurfaceValue) {
2248 float temp = SurfaceValue[v1];
2249 temp /= MaxSurfaceValue;
2250 SetSurfaceColor(temp, SurfaceColor, NColor, red, green, blue);
2251 glColor4f(red, green, blue, alpha);
2252 }
2253
2254 glVertex3d(Vertices[v1].x, Vertices[v1].y, Vertices[v1].z);
2255
2256 if (SurfaceValue) {
2257 float temp = SurfaceValue[v2];
2258 temp /= MaxSurfaceValue;
2259 SetSurfaceColor(temp, SurfaceColor, NColor, red, green, blue);
2260 glColor4f(red, green, blue, alpha);
2261 }
2262
2263 glVertex3d(Vertices[v2].x, Vertices[v2].y, Vertices[v2].z);
2264 //Once for 1-2, once for 2-3
2265 glVertex3d(Vertices[v2].x, Vertices[v2].y, Vertices[v2].z);
2266
2267 if (SurfaceValue) {
2268 float temp = SurfaceValue[v3];
2269 temp /= MaxSurfaceValue;
2270 SetSurfaceColor(temp, SurfaceColor, NColor, red, green, blue);
2271 glColor4f(red, green, blue, alpha);
2272 }
2273 glVertex3d(Vertices[v3].x, Vertices[v3].y, Vertices[v3].z);
2274
2275 glVertex3d(Vertices[v3].x, Vertices[v3].y, Vertices[v3].z);
2276
2277 if (SurfaceValue) {
2278 float temp = SurfaceValue[v1];
2279 temp /= MaxSurfaceValue;
2280 SetSurfaceColor(temp, SurfaceColor, NColor, red, green, blue);
2281 glColor4f(red, green, blue, alpha);
2282 }
2283
2284 glVertex3d(Vertices[v1].x, Vertices[v1].y, Vertices[v1].z);
2285 }
2286 glEnd(); //End of Line creation
2287 return result;
2288 }
2289
2290 //Utility function to create a line made up of a variable width cylinder. the GLUquadricObj must be preallocated
CreateCylinderFromLine(GLUquadricObj * qobj,const CPoint3D & lineStart,const CPoint3D & lineEnd,float lineWidth,int nslices,int nstacks,bool cap)2291 void CreateCylinderFromLine(GLUquadricObj *qobj, const CPoint3D& lineStart, const CPoint3D& lineEnd,
2292 float lineWidth, int nslices, int nstacks, bool cap) {
2293
2294 if (qobj == NULL) return;
2295
2296 CPoint3D offset;
2297 CPoint3D NormalOffset(0.0f, 0.0f, 0.0f);
2298 CPoint3D NormStart(0.0f, 0.0f, 1.0f);
2299 Matrix4D rotMat;
2300
2301 offset = lineEnd - lineStart;
2302 float length = offset.Magnitude();
2303 if (length>0.00001) {
2304 NormalOffset.x = offset.x/length;
2305 NormalOffset.y = offset.y/length;
2306 NormalOffset.z = offset.z/length;
2307 }
2308
2309 SetRotationMatrix(rotMat, &NormStart, &NormalOffset);
2310 rotMat[3][0] = lineStart.x;
2311 rotMat[3][1] = lineStart.y;
2312 rotMat[3][2] = lineStart.z;
2313
2314 glPushMatrix();
2315 glMultMatrixf((const GLfloat *) &rotMat);
2316 gluCylinder(qobj, lineWidth, lineWidth, length, nslices, nstacks);
2317 if (cap) {
2318 glTranslatef(0.0, 0.0, length);
2319 gluSphere(qobj, lineWidth, nslices, nstacks);
2320 /* gluDisk(qobj, 0.0, lineWidth, nslices, 2); */
2321 }
2322 glPopMatrix();
2323 }
2324
DashedQuadFromLine(const CPoint3D & pt1,const CPoint3D & pt2,float width,float m[16],const CPoint3D &,float offset,GLuint length_anno_tex_id,const WinPrefs * Prefs,bool draw_label)2325 void DashedQuadFromLine(const CPoint3D& pt1,
2326 const CPoint3D& pt2, float width, float m[16],
2327 const CPoint3D& /*x_world*/, float offset,
2328 GLuint length_anno_tex_id,
2329 const WinPrefs * Prefs, bool draw_label) {
2330
2331 float len;
2332 GLdouble scr_coords1[3]; // Screen coordinates of pt1
2333 GLdouble scr_coords2[3]; // Screen coordinates of pt2
2334 CPoint3D scr_vec; // Screen space vector between atoms
2335 GLdouble perp_obj[3]; // Object coords on vector perp. to scr_vec
2336 CPoint3D offset_vec; // Direction to shift bond cylinders
2337 GLdouble modelview[16];
2338 GLdouble proj[16];
2339 GLint viewport[4];
2340
2341 glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
2342 glGetDoublev(GL_PROJECTION_MATRIX, proj);
2343 glGetIntegerv(GL_VIEWPORT, viewport);
2344
2345 // Find screen coordinates of one point.
2346 gluProject(pt1.x, pt1.y, pt1.z, modelview, proj, viewport,
2347 &(scr_coords1[0]), &(scr_coords1[1]), &(scr_coords1[2]));
2348
2349 // Find screen coordinates of other atom.
2350 gluProject(pt2.x, pt2.y, pt2.z, modelview, proj, viewport,
2351 &(scr_coords2[0]), &(scr_coords2[1]), &(scr_coords2[2]));
2352
2353 // Find vector perpendicular to vector between two screen points and
2354 // normalize it so we can scalar multiply it later. We flip and
2355 // negate the slope of the line between the two screen coordinates to
2356 // get the slope of the perpendicular line.
2357 scr_vec.x = scr_coords2[1] - scr_coords1[1];
2358 scr_vec.y = scr_coords1[0] - scr_coords2[0];
2359 scr_vec.z = 0;
2360 scr_vec *= 1 / scr_vec.Magnitude();
2361
2362 // Now find a point on the perpendicular vector with pt1's depth
2363 // and get its object coordinates.
2364 gluUnProject(scr_coords1[0] + scr_vec.x * 10,
2365 scr_coords1[1] + scr_vec.y * 10,
2366 scr_coords1[2],
2367 modelview, proj, viewport,
2368 &(perp_obj[0]), &(perp_obj[1]), &(perp_obj[2]));
2369
2370 // Finally, we see what direction all bond cylinders must be offset
2371 // so that they will always stay in view.
2372 offset_vec.x = perp_obj[0] - pt1.x;
2373 offset_vec.y = perp_obj[1] - pt1.y;
2374 offset_vec.z = perp_obj[2] - pt1.z;
2375 offset_vec *= 1 / offset_vec.Magnitude();
2376
2377 len = (pt2 - pt1).Magnitude();
2378
2379 CPoint3D new_pt1a;
2380 CPoint3D new_pt1b;
2381 CPoint3D new_pt2a;
2382 CPoint3D new_pt2b;
2383
2384 new_pt1a.x = pt1.x + offset_vec.x * width;
2385 new_pt1a.y = pt1.y + offset_vec.y * width;
2386 new_pt1a.z = pt1.z + offset_vec.z * width;
2387 new_pt1b.x = pt1.x - offset_vec.x * width;
2388 new_pt1b.y = pt1.y - offset_vec.y * width;
2389 new_pt1b.z = pt1.z - offset_vec.z * width;
2390 new_pt2a.x = pt2.x + offset_vec.x * width;
2391 new_pt2a.y = pt2.y + offset_vec.y * width;
2392 new_pt2a.z = pt2.z + offset_vec.z * width;
2393 new_pt2b.x = pt2.x - offset_vec.x * width;
2394 new_pt2b.y = pt2.y - offset_vec.y * width;
2395 new_pt2b.z = pt2.z - offset_vec.z * width;
2396
2397 glEnable(GL_ALPHA_TEST);
2398 glAlphaFunc(GL_GREATER, 0.5f);
2399
2400 glBindTexture(GL_TEXTURE_2D, length_anno_tex_id);
2401 glEnable(GL_TEXTURE_2D);
2402 glBegin(GL_QUADS);
2403 glTexCoord2f(0.0f, 1.0f);
2404 glVertex3f(new_pt1a.x, new_pt1a.y, new_pt1a.z);
2405 glTexCoord2f(0.0f, 0.0f);
2406 glVertex3f(new_pt1b.x, new_pt1b.y, new_pt1b.z);
2407 glTexCoord2f(len / (5.0f * width), 0.0f);
2408 glVertex3f(new_pt2b.x, new_pt2b.y, new_pt2b.z);
2409 glTexCoord2f(len / (5.0f * width), 1.0f);
2410 glVertex3f(new_pt2a.x, new_pt2a.y, new_pt2a.z);
2411 glEnd();
2412 glDisable(GL_TEXTURE_2D);
2413 glDisable(GL_ALPHA_TEST);
2414
2415 if (draw_label) {
2416 char len_label[40];
2417 float LabelSize = Prefs->GetAnnotationLabelSize();
2418
2419 // We move the midpoint of the line and align the viewer to draw the
2420 // text.
2421 glPushMatrix();
2422 glTranslatef((pt1.x + pt2.x) / 2.0f,
2423 (pt1.y + pt2.y) / 2.0f,
2424 (pt1.z + pt2.z) / 2.0f);
2425 glMultMatrixf(m);
2426
2427 // Move out for the kind of bond that exists between the atoms.
2428 glTranslatef(-offset - 0.01f, 0.0f, 0.0f);
2429
2430 glScalef(-5.0f * LabelSize, 5.0f * LabelSize, 0.1f);
2431
2432 // This apparently is some magic number. GLF doesn't start drawing the
2433 // string at the origin, doesn't consider the anchor point, and doesn't
2434 // consider the bounding box it returns. So, we have a magic number.
2435 // It is indirectly dependent on LabelSize from the scaling above, so
2436 // it looks constant even though it is not.
2437 glTranslatef(0.004f, 0.0f, 0.0f);
2438
2439 sprintf(len_label, "%.2f", len);
2440 DrawString(len_label);
2441
2442 glPopMatrix();
2443 }
2444
2445 }
2446
DrawRotationAxis(const CPoint3D & lineStart,const CPoint3D & lineEnd,const int & order)2447 void DrawRotationAxis(const CPoint3D & lineStart, const CPoint3D & lineEnd, const int & order) {
2448
2449 float plane_emissive[] = { 0.0, 0.3, 0.7, 0.2 };
2450 // float plane_diffuse[] = { 0.0, 0.3, 0.6, 0.3 };
2451 float plane_specular[] = { 0.0, 0.3, 0.6, 1.0 };
2452
2453 glPushAttrib(GL_LIGHTING_BIT);
2454
2455 int imageWidth =16;
2456 //Our width needs to be a power of two. So orders 1, 2 and 4 are no problem. Other orders such as 3
2457 //must be padded out to the next higher power of two.
2458 int repeat = order;
2459 if (order == 3) repeat = 4;
2460 if (order > 4) repeat = 8;
2461 GLubyte bw[16][16] ={ {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
2462 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
2463 {0,0,0,0,0,0,0,0,255,0,0,0,0,0,0,0},
2464 {0,0,0,0,0,0,0,0,255,128,0,0,0,0,0,0},
2465 {0,0,0,0,0,0,0,0,255,255,0,0,0,0,0,0},
2466 {0,255,255,255,255,255,255,255,255,255,128,0,0,0,0,0},
2467 {0,255,255,255,255,255,255,255,255,255,255,0,0,0,0},
2468 {0,255,255,255,255,255,255,255,255,255,255,128,0,0,0},
2469 {0,255,255,255,255,255,255,255,255,255,255,128,0,0,0},
2470 {0,255,255,255,255,255,255,255,255,255,255,0,0,0,0},
2471 {0,255,255,255,255,255,255,255,255,255,128,0,0,0,0,0},
2472 {0,0,0,0,0,0,0,0,255,255,0,0,0,0,0,0},
2473 {0,0,0,0,0,0,0,0,255,128,0,0,0,0,0,0},
2474 {0,0,0,0,0,0,0,0,255,0,0,0,0,0,0,0},
2475 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
2476 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}};
2477 GLubyte * testimage = new GLubyte[imageWidth*imageWidth*repeat*4];
2478 int padding = ((repeat - order)*16)/order;
2479 int remainder = ((repeat - order)*16)-(padding*order);
2480 int p = 0;
2481 for (int i=0; i<imageWidth; i++) {
2482 for (int o=0; o<order; o++) {
2483 for (int j=0; j<imageWidth; j++) {
2484 testimage[p] = bw[i][j];
2485 testimage[p + 1] = bw[i][j];
2486 testimage[p + 2] = bw[i][j];
2487 testimage[p + 3] = bw[i][j];
2488 if (bw[i][j] == 0) testimage[p + 3]=178;
2489 p+=4;
2490 }
2491 for (int t=0; t<padding; t++) {
2492 testimage[p] = 0;
2493 testimage[p + 1] = 0;
2494 testimage[p + 2] = 0;
2495 testimage[p + 3] = 178;
2496 p+=4;
2497 }
2498 }
2499 for (int t=0; t<remainder; t++) {
2500 testimage[p] = 0;
2501 testimage[p + 1] = 0;
2502 testimage[p + 2] = 0;
2503 testimage[p + 3] = 178;
2504 p+=4;
2505 }
2506 }
2507 glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
2508 GLuint texname;
2509 glGenTextures(1, &texname);
2510 glBindTexture(GL_TEXTURE_2D, texname);
2511
2512 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
2513 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
2514 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
2515 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
2516
2517 glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, imageWidth*repeat, imageWidth, 0, GL_RGBA,
2518 GL_UNSIGNED_BYTE, testimage);
2519
2520 glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, plane_emissive);
2521 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, plane_specular);
2522 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 30.0f);
2523 glColor4f(0, .64, .85, 0.7);
2524 /* glColor4fv(plane_diffuse); */
2525
2526 glEnable(GL_BLEND);
2527 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
2528 glEnable(GL_TEXTURE_2D);
2529 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_BLEND);
2530
2531 GLUquadricObj * qobj = NULL;
2532 qobj = gluNewQuadric();
2533
2534 CPoint3D offset, NormalOffset, NormEnd, NormStart = CPoint3D(0.0f, 0.0f, 1.0f);
2535 Matrix4D rotMat;
2536
2537 offset.x = lineEnd.x - lineStart.x;
2538 offset.y = lineEnd.y - lineStart.y;
2539 offset.z = lineEnd.z - lineStart.z;
2540 float length = offset.Magnitude();
2541 if (length>0.00001) {
2542 NormalOffset.x = offset.x/length;
2543 NormalOffset.y = offset.y/length;
2544 NormalOffset.z = offset.z/length;
2545 } else
2546 NormalOffset.x=NormalOffset.y=NormalOffset.z=0.0;
2547 NormEnd = lineEnd;
2548 Normalize3D(&NormEnd);
2549 SetRotationMatrix(rotMat, &NormStart, &NormalOffset);
2550 rotMat[3][0] = lineStart.x;
2551 rotMat[3][1] = lineStart.y;
2552 rotMat[3][2] = lineStart.z;
2553
2554 glPushMatrix();
2555 glMultMatrixf((const GLfloat *) &rotMat);
2556 gluQuadricTexture(qobj, GL_TRUE);
2557 gluCylinder(qobj, 0.1, 0.1, 0.2, 12, 1);
2558 glTranslatef(0, 0, 0.2);
2559 glDisable(GL_TEXTURE_2D);
2560 gluQuadricTexture(qobj, GL_FALSE);
2561 gluCylinder(qobj, 0.1, 0.1, length-0.4, 12, 1);
2562 glEnable(GL_TEXTURE_2D);
2563 glTranslatef(0, 0, length-0.4);
2564 gluQuadricTexture(qobj, GL_TRUE);
2565 gluCylinder(qobj, 0.1, 0.1, 0.2, 12, 1);
2566 glPopMatrix();
2567
2568 if (qobj) gluDeleteQuadric(qobj); //finally delete the quadric object
2569
2570 glDisable(GL_TEXTURE_2D);
2571 glDisable(GL_BLEND);
2572 glDeleteTextures(1, &texname);
2573
2574 glPopAttrib();
2575
2576 delete[] testimage;
2577 }
2578
2579 //Draw a translucent plane
DrawTranslucentPlane(const CPoint3D & origin,const CPoint3D & p1,const CPoint3D & p2)2580 void DrawTranslucentPlane(const CPoint3D & origin, const CPoint3D & p1, const CPoint3D & p2) {
2581
2582 float plane_emissive[] = { 0.0, 0.3, 0.7, 0.2 };
2583 // float plane_diffuse[] = { 0.0, 0.3, 0.6, 0.3 };
2584 float plane_specular[] = { 0.0, 0.3, 0.6, 1.0 };
2585 glPushAttrib(GL_LIGHTING_BIT);
2586 Matrix4D rotationMatrix;
2587
2588 CPoint3D vec1 = p1 - origin;
2589 CPoint3D vec2 = p2 - origin;
2590 float s1Length = vec1.Magnitude();
2591 float s2Length = vec2.Magnitude();
2592 SetPlaneRotation(rotationMatrix, vec1, vec2);
2593
2594 glPushMatrix();
2595 glTranslatef(origin.x, origin.y, origin.z);
2596 glMultMatrixf((const GLfloat *) &rotationMatrix);
2597
2598 // Symmetry planes are viewed from either side, so if we should either
2599 // light both sides or light neither. Neither is simpler.
2600 glDisable(GL_LIGHTING);
2601
2602 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
2603 glEnable(GL_BLEND);
2604 glDepthMask(GL_FALSE);
2605
2606 glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, plane_emissive);
2607 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, plane_specular);
2608 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 30.0f);
2609 glColor4f(0, .64, .85, 0.3);
2610 glRectf(0, 0, s1Length, s2Length);
2611 glDisable(GL_BLEND);
2612
2613 glColor4f(0, .64, .85, 1);
2614 glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
2615 glLineWidth(1);
2616 glRectf(0, 0, s1Length, s2Length);
2617 glEnable(GL_LIGHTING);
2618 glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
2619
2620 glDepthMask(GL_TRUE);
2621 glPopAttrib();
2622
2623 glPopMatrix();
2624 }
2625
DrawInversionPoint(void)2626 void DrawInversionPoint(void) {
2627 GLUquadricObj * qobj = NULL;
2628 qobj = gluNewQuadric();
2629 // float sphere_emissive[] = { 0.0, 0.3, 0.7, 0.2 };
2630 // float sphere_diffuse[] = { 1.0, 1.0, 1.0, 1.0 };
2631 float sphere_specular[] = { 1.0, 1.0, 1.0, 1.0 };
2632
2633 glPushAttrib(GL_LIGHTING_BIT);
2634
2635 glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 128.0);
2636 // glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, sphere_emissive);
2637 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, sphere_specular);
2638 glColor4f(.7, .7, .7, 1.0);
2639
2640 //Assume the inversion point is at the origin
2641 gluSphere(qobj, 0.2, 30, 20); //Create and draw the sphere
2642
2643 //Add several arrows to indicate the inversion...
2644 float plane_emissive[] = { 0.0, 0.3, 0.7, 0.2 };
2645 // float plane_diffuse[] = { 0.0, 0.3, 0.6, 0.3 };
2646 float plane_specular[] = { 0.0, 0.3, 0.6, 1.0 };
2647 glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, plane_emissive);
2648 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, plane_specular);
2649 glColor4f(0, .64, .85, 0.7);
2650
2651 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
2652 glEnable(GL_BLEND);
2653
2654 gluCylinder(qobj, 0.02, 0.02, 0.6, 12, 1);
2655 glPushMatrix();
2656 glTranslatef(0.0, 0.0, 0.6);
2657 gluDisk(qobj, 0.0, 0.06, 12, 2);
2658 gluCylinder(qobj, 0.06, 0.0, 0.07, 12, 3);
2659 glPopMatrix();
2660
2661 glPushMatrix();
2662 Matrix4D rotMat;
2663 CPoint3D vector = CPoint3D(1.0f, 0.0f, 0.0), NormStart = CPoint3D(0.0f, 0.0f,1.0f);
2664 InitRotationMatrix(rotMat);
2665 rotMat[2][2] = -1;
2666 glMultMatrixf((const GLfloat *) &rotMat);
2667 DrawArrow(0.7, 0.02, 12);
2668 glPopMatrix();
2669 glPushMatrix();
2670 SetRotationMatrix(rotMat, &NormStart, &vector);
2671 glMultMatrixf((const GLfloat *) &rotMat);
2672 DrawArrow(0.7, 0.02, 12);
2673 glPopMatrix();
2674 glPushMatrix();
2675 rotMat[2][0] *= -1.0;
2676 glMultMatrixf((const GLfloat *) &rotMat);
2677 DrawArrow(0.7, 0.02, 12);
2678 glPopMatrix();
2679 vector.x = 0.0;
2680 vector.y = 1.0;
2681 glPushMatrix();
2682 SetRotationMatrix(rotMat, &NormStart, &vector);
2683 glMultMatrixf((const GLfloat *) &rotMat);
2684 DrawArrow(0.7, 0.02, 12);
2685 glPopMatrix();
2686 glPushMatrix();
2687 rotMat[2][1] *= -1.0;
2688 glMultMatrixf((const GLfloat *) &rotMat);
2689 DrawArrow(0.7, 0.02, 12);
2690 glPopMatrix();
2691
2692 glDisable(GL_BLEND);
2693 glPopAttrib();
2694 if (qobj) gluDeleteQuadric(qobj); //finally delete the quadric object
2695 }
2696
DrawArrow(const float & length,const float & width,const int & quality)2697 void DrawArrow(const float & length, const float & width, const int & quality) {
2698 //Draw an arrow with arrow head at far end.
2699
2700 GLUquadricObj * qobj = NULL;
2701 qobj = gluNewQuadric();
2702 gluDisk(qobj, 0.0, width, quality, 2);
2703 gluCylinder(qobj, width, width, length-3.5*width, quality, 1);
2704 glPushMatrix();
2705 glTranslatef(0.0, 0.0, length-2.5*width);
2706 gluDisk(qobj, 0.0, 3*width, quality, 2);
2707 gluCylinder(qobj, 3*width, 0.0, 3.5*width, quality, 2);
2708 glPopMatrix();
2709 if (qobj) gluDeleteQuadric(qobj); //finally delete the quadric object
2710 }
2711
2712 /**
2713 * This draws a dashed angle between the three specified points and
2714 * marks the angle's degrees with a text label. The angle goes from
2715 * pt1 to pt2 to pt3. To complicate things a bit, we draw each dash
2716 * so that it is aligned to the screen so the viewer can always see it.
2717 */
DrawAngleAnnotation(const CPoint3D * pt1,const CPoint3D * pt2,const CPoint3D * pt3,const WinPrefs * Prefs,GLuint length_anno_tex_id,CPoint3D * ambig_axis)2718 void DrawAngleAnnotation(const CPoint3D *pt1, const CPoint3D *pt2,
2719 const CPoint3D *pt3, const WinPrefs *Prefs,
2720 GLuint length_anno_tex_id, CPoint3D *ambig_axis) {
2721
2722 CPoint3D vec1;
2723 CPoint3D vec2;
2724 Matrix4D plane2xy;
2725 float min_len;
2726 float len1;
2727 float len2;
2728 float angle;
2729 float m[16];
2730 float chord_len;
2731 float x_eye[3];
2732 float y_eye[3];
2733 float delta;
2734
2735 // First we get vectors representing the angle.
2736 vec1 = *pt1 - *pt2;
2737 vec2 = *pt3 - *pt2;
2738
2739 // We need to find the shorter of the two vectors so that we know how far
2740 // along each vector we can draw the annotation. And with the lengths
2741 // saved, we can go ahead and normalize.
2742 len1 = vec1.Magnitude();
2743 len2 = vec2.Magnitude();
2744 min_len = MIN(len1, len2);
2745
2746 vec1 = vec1 * (1.0f / len1);
2747 vec2 = vec2 * (1.0f / len2);
2748
2749 // Grab the angle.
2750 float dot = DotProduct3D(&vec1, &vec2);
2751 angle = acos(dot);
2752
2753 // If two vectors are not (anti-)parallel, drawing the annotation is easy.
2754 // We move into the space where vec1 is the x-axis and vec2 is in the
2755 // xy-plane. In this space we can draw a circular segment with simple
2756 // geometry.
2757 if (fabs(dot) < 0.999999f) {
2758 SetPlaneRotation(plane2xy, vec1, vec2);
2759 }
2760
2761 // If they are (anti-)parallel, there're many possible planes we can move
2762 // into.
2763 else {
2764
2765 // acos doesn't tolerate values whose absolute values are > 1, which
2766 // might happen because of precision errors. So, we fix the angle.
2767 angle = dot > 0.0f ? 0.0f : kPi;
2768
2769 // For dihedral annotations, we push the normal of the dihedral planes
2770 // into the new space's xy-plane.
2771 if (ambig_axis) {
2772 SetPlaneRotation(plane2xy, vec1, *ambig_axis);
2773 }
2774
2775 // For angle annotations, we try to keep them as visible as possible.
2776 // We keep vec1 as the x-axis in our new space, and we try to put
2777 // either the screen's x- or y-axis in the new space's xy-plane.
2778 else {
2779 glGetFloatv(GL_MODELVIEW_MATRIX, m);
2780 CPoint3D vec3;
2781 vec3.x = m[0];
2782 vec3.y = m[4];
2783 vec3.z = m[8];
2784
2785 if (fabs(DotProduct3D(&vec1, &vec2)) < 0.999999f) {
2786 SetPlaneRotation(plane2xy, vec1, vec3);
2787 } else {
2788 vec3.x = m[1];
2789 vec3.y = m[5];
2790 vec3.z = m[9];
2791 SetPlaneRotation(plane2xy, vec1, vec3);
2792 }
2793 }
2794 }
2795
2796 glPushMatrix();
2797
2798 // We move to the center vertex and transform into the space described
2799 // above.
2800 glTranslatef(pt2->x, pt2->y, pt2->z);
2801 glMultMatrixf((const GLfloat *) &plane2xy);
2802
2803 // We want screen-aligned dashes. So we want to draw quads that extend in
2804 // the x- and y-directions in eye space. Since we draw in object space, we
2805 // must figure out what the eye's x- and y- axes are there. To do that,
2806 // all we have to do is grab the first two axes from the inverse of the
2807 // modelview matrix.
2808 glGetFloatv(GL_MODELVIEW_MATRIX, m);
2809
2810 x_eye[0] = m[0];
2811 x_eye[1] = m[4];
2812 x_eye[2] = m[8];
2813
2814 y_eye[0] = m[1];
2815 y_eye[1] = m[5];
2816 y_eye[2] = m[9];
2817
2818 GLdouble BondSize = Prefs->GetQD3DBondWidth();
2819 CPoint3D origin(0.0f, 0.0f, 0.0f);
2820
2821 // How long should each dash be, and how often do we want them?
2822 chord_len = BondSize * 0.25f * 10.0f;
2823 delta = angle / (min_len * angle / chord_len);
2824
2825 float ca = cos(angle);
2826 float sa = sin(angle);
2827
2828 // First we draw the dashes along each side.
2829 DashedQuadFromLine(origin, CPoint3D(min_len, 0.0f, 0.0f),
2830 BondSize * 0.25f, m,
2831 CPoint3D(x_eye[0], x_eye[1], x_eye[2]),
2832 0.0f, length_anno_tex_id, Prefs, false);
2833 DashedQuadFromLine(CPoint3D(min_len * ca, min_len * sa, 0.0f), origin,
2834 BondSize * 0.25f, m,
2835 CPoint3D(x_eye[0], x_eye[1], x_eye[2]),
2836 0.0f, length_anno_tex_id, Prefs, false);
2837
2838 // Then we draw a series of chords that approximates the angle's arc.
2839 float next_delta;
2840 for (float i = 0.0f; i < angle; i += delta) {
2841
2842 // We want the chord to go from the current slice of the pie to the
2843 // next, but since the slice angle is constant and the anno angle may
2844 // not be a multiple of the slice angle, we don't want to overshoot the
2845 // anno angle. Also, we subtract off an epsilon so that the texture
2846 // doesn't start to repeat.
2847 next_delta = MIN(i + delta, angle) - 0.01f;
2848
2849 DashedQuadFromLine(CPoint3D(min_len * cos(i), min_len * sin(i), 0.0f),
2850 CPoint3D(min_len * cos(next_delta), min_len * sin(next_delta), 0.0f),
2851 BondSize * 0.25f, m,
2852 CPoint3D(x_eye[0], x_eye[1], x_eye[2]),
2853 0.0f, length_anno_tex_id, Prefs, false);
2854 }
2855
2856 CPoint3D lookat_eye = CPoint3D(0.0f, 0.0f, 1.0f);
2857 CPoint3D up_eye = CPoint3D(0.0f, 1.0f, 0.0f);
2858 CPoint3D lookat_world;
2859 CPoint3D up_world;
2860 CPoint3D r;
2861
2862 Matrix4D mv_inv;
2863
2864 // Invert just the rotation portion of the modelview matrix. We can't
2865 // use InverseMatrix because it considers the translation factors.
2866 // We don't want the translation factors because we want to transform
2867 // a vector, and Rotate3DPt transforms only points (and not vectors).
2868 mv_inv[0][0] = m[0];
2869 mv_inv[0][1] = m[4];
2870 mv_inv[0][2] = m[8];
2871 mv_inv[0][3] = 0.0f;
2872 mv_inv[1][0] = m[1];
2873 mv_inv[1][1] = m[5];
2874 mv_inv[1][2] = m[9];
2875 mv_inv[1][3] = 0.0f;
2876 mv_inv[2][0] = m[2];
2877 mv_inv[2][1] = m[6];
2878 mv_inv[2][2] = m[10];
2879 mv_inv[2][3] = 0.0f;
2880 mv_inv[3][0] = 0.0f;
2881 mv_inv[3][1] = 0.0f;
2882 mv_inv[3][2] = 0.0f;
2883 mv_inv[3][3] = 1.0f;
2884
2885 // Transform the eye space vectors to world coordinates, and find
2886 // a third vector to form a basis set.
2887 Rotate3DPt(mv_inv, lookat_eye, &lookat_world);
2888 Rotate3DPt(mv_inv, up_eye, &up_world);
2889 CrossProduct3D(&lookat_world, &up_world, &r);
2890
2891 m[0] = r.x;
2892 m[1] = r.y;
2893 m[2] = r.z;
2894 m[3] = 0.0f;
2895
2896 m[4] = up_world.x;
2897 m[5] = up_world.y;
2898 m[6] = up_world.z;
2899 m[7] = 0.0f;
2900
2901 m[8] = lookat_world.x;
2902 m[9] = lookat_world.y;
2903 m[10] = lookat_world.z;
2904 m[11] = 0.0f;
2905
2906 m[12] = m[13] = m[14] = 0.0f;
2907 m[15] = 1.0f;
2908
2909 char angle_label[40];
2910 sprintf(angle_label, "%.2f", angle * 180.f / kPi);
2911 glTranslatef(min_len * cos(angle * 0.5f),
2912 min_len * sin(angle * 0.5f), 0.0f);
2913 glMultMatrixf(m);
2914 float LabelSize = Prefs->GetAnnotationLabelSize();
2915 glTranslatef(-0.1f + chord_len, 0.0f, 0.0f);
2916 glScalef(-5.0f * LabelSize, 5.0f * LabelSize, 0.1f);
2917
2918 // glfDrawSolidString(angle_label);
2919 DrawString(angle_label);
2920
2921 glPopMatrix();
2922
2923 }
2924
DrawSceneString(const float scale_factor,const float shift_x,const float shift_y,const float shift_z,const wxString & label)2925 void DrawSceneString(const float scale_factor, const float shift_x,
2926 const float shift_y, const float shift_z,
2927 const wxString& label) {
2928
2929 // This function draws a string (label) at the origin of the current
2930 // coordinate system, offset by shift_[xyz] and aligned with the screen.
2931 // Right before the string is drawn, the coordinate system is uniformly
2932 // scaled by scale_factor (since glf draws in a larger coordinate system
2933 // than we want).
2934 //
2935 // It's assumed that the color and coordinate system are in the correct
2936 // state when this function is called. It's also assumed that lighting is
2937 // turned on. The assumed state is restored before the function returns.
2938
2939 float m[16]; // OpenGL matrix
2940 Matrix4D mv_inv; // Inverse of modelview matrix
2941 CPoint3D lookat_eye = CPoint3D(0.0f, 0.0f, 1.0f); // Lookat vector in eye-space
2942 CPoint3D lookat_world; // transformed to world-space
2943 CPoint3D up_eye = CPoint3D(0.0f, 1.0f, 0.0f); // Up vector in eye-space
2944 CPoint3D up_world; // transformed to world-space
2945 CPoint3D r; // Right vector in world-space
2946
2947 // Invert just the rotation portion of the modelview matrix. This is
2948 // much faster than inverting an arbitrary matrix.
2949 glGetFloatv(GL_MODELVIEW_MATRIX, m);
2950 mv_inv[0][0] = m[0];
2951 mv_inv[0][1] = m[4];
2952 mv_inv[0][2] = m[8];
2953 mv_inv[0][3] = 0.0f;
2954 mv_inv[1][0] = m[1];
2955 mv_inv[1][1] = m[5];
2956 mv_inv[1][2] = m[9];
2957 mv_inv[1][3] = 0.0f;
2958 mv_inv[2][0] = m[2];
2959 mv_inv[2][1] = m[6];
2960 mv_inv[2][2] = m[10];
2961 mv_inv[2][3] = 0.0f;
2962 mv_inv[3][0] = 0.0f;
2963 mv_inv[3][1] = 0.0f;
2964 mv_inv[3][2] = 0.0f;
2965 mv_inv[3][3] = 1.0f;
2966
2967 // Transform the eye space vectors to world coordinates, and find
2968 // a third vector to form a basis set.
2969 Rotate3DPt(mv_inv, lookat_eye, &lookat_world);
2970 Rotate3DPt(mv_inv, up_eye, &up_world);
2971 CrossProduct3D(&lookat_world, &up_world, &r);
2972
2973 // Enter the basis into a new modelview matrix.
2974 m[0] = r.x;
2975 m[1] = r.y;
2976 m[2] = r.z;
2977 m[3] = 0.0f;
2978
2979 m[4] = up_world.x;
2980 m[5] = up_world.y;
2981 m[6] = up_world.z;
2982 m[7] = 0.0f;
2983
2984 m[8] = lookat_world.x;
2985 m[9] = lookat_world.y;
2986 m[10] = lookat_world.z;
2987 m[11] = 0.0f;
2988
2989 m[12] = m[13] = m[14] = 0.0f;
2990 m[15] = 1.0f;
2991
2992 glDisable(GL_LIGHTING);
2993
2994 // Undo the rotation by multiplying by the inverse. Then shift, scale, and
2995 // draw the string.
2996 glPushMatrix();
2997 glMultMatrixf(m);
2998 glTranslatef(shift_x, shift_y, shift_z);
2999 glScalef(-scale_factor, scale_factor, scale_factor);
3000 glfDrawSolidString(label.mb_str());
3001 glPopMatrix();
3002
3003 glEnable(GL_LIGHTING);
3004
3005 }
3006
DrawString(const char * str)3007 void DrawString(const char *str) {
3008 // glfDrawSolidString(str);
3009
3010 glfStartBitmapDrawing();
3011
3012 glBlendFunc(GL_SRC_ALPHA, GL_DST_ALPHA);
3013 glEnable(GL_BLEND);
3014 glfDrawBMaskString(str);
3015 glDisable(GL_BLEND);
3016
3017 glfStopBitmapDrawing();
3018 }
3019
3020 #define CYL_RADIUS 0.05f
3021 #include <algorithm>
DrawBondingSites(long iatom,float radius,GLUquadricObj * qobj,int site_id,CPoint3D * vector)3022 void MolDisplayWin::DrawBondingSites(long iatom, float radius, GLUquadricObj *qobj, int site_id, CPoint3D * vector) {
3023
3024 // This function either draws remaining bonding sites for the specified
3025 // atom or returns the vector along the atom's specified bonding site. If
3026 // site_vec is 0, drawing is done. Otherwise, the normalized vector for
3027 // the site_vec is returned.
3028
3029 Frame *lFrame=MainData->cFrame;
3030 mpAtom *lAtoms = lFrame->Atoms;
3031 Bond *lBonds = lFrame->Bonds;
3032 long NumBonds = lFrame->NumBonds;
3033
3034 short coordination = lAtoms[iatom].GetCoordinationNumber();
3035 std::vector<Bond *> bonds;
3036 int bonded_atoms[6];
3037 CPoint3D vecs[6];
3038 unsigned int i;
3039
3040 bonds.reserve(6);
3041
3042 if (coordination == 0) return;
3043
3044 // First we grab all bonds pertaining to the atom. If the number of bonds
3045 // exceeds five, then there's no point in showing any bonding sites since
3046 // the atom is either full or overfull.
3047 for (i = 0; i < NumBonds; i++) {
3048 if (((iatom == lBonds[i].Atom1) || (iatom == lBonds[i].Atom2)) &&
3049 (lBonds[i].Order > kHydrogenBond)) {
3050 bonds.push_back(&lBonds[i]);
3051 }
3052 }
3053
3054 if (bonds.size() >= coordination) {
3055 return;
3056 }
3057
3058 // Test to see if any bond sites are unused, if not return.
3059 // if (bonds.size() >= coordination) return;
3060
3061 // We sort the bonds by bond order, and figure out what the atom opposite
3062 // this one is. We'll need this info to find the vector between the two
3063 // atoms.
3064 std::sort(bonds.begin(), bonds.end(), Bond::RefLessThan);
3065 for (i = 0; i < bonds.size(); i++) {
3066 bonded_atoms[i] = (iatom == bonds[i]->Atom1) ?
3067 bonds[i]->Atom2 : bonds[i]->Atom1;
3068 }
3069
3070 // Calculate vectors from central atom to each bonded atom.
3071 for (i = 0; i < bonds.size(); i++) {
3072 vecs[i] = lAtoms[bonded_atoms[i]].Position - lAtoms[iatom].Position;
3073 Normalize3D(&vecs[i]);
3074 }
3075
3076 int lpCount = lAtoms[iatom].GetLonePairCount();
3077
3078 CPoint3D origin = CPoint3D(0.0f, 0.0f, 0.0f);
3079 float c, s, b;
3080
3081 Matrix4D rot;
3082 InitRotationMatrix(rot);
3083 CPoint3D b1Offset(0.0f, 1.0f, 0.0f);
3084 if (bonds.size() > 0) { //orient our frame to the primary bond
3085 b1Offset = lAtoms[bonded_atoms[0]].Position - lAtoms[iatom].Position;
3086 Normalize3D(&b1Offset);
3087 CPoint3D target(0.0f, 1.0f, 0.0f);
3088 SetRotationMatrix(rot, &target, &b1Offset);
3089 }
3090
3091 if (site_id == 0) {
3092 glPushName(0);
3093 glColor3f(0.9f, 0.1f, 0.1f);
3094 }
3095
3096 #define DO_SITE(vec, vec_site) \
3097 if (site_id == 0) { \
3098 glLoadName(vec_site); \
3099 CreateCylinderFromLine(qobj, origin, \
3100 vec * 2.0f * radius, \
3101 CYL_RADIUS, 10, 3, true); \
3102 } else if (site_id == vec_site) { \
3103 *vector = vec; \
3104 return; \
3105 }
3106
3107 // If we get here, we know that there are some unbonded sites.
3108 switch (coordination+lpCount) {
3109
3110 // Configuration will be linear.
3111 case 1:
3112 {
3113 CPoint3D v(0.0f, 1.0f, 0.0f);
3114 DO_SITE(v, 1);
3115 }
3116 break;
3117
3118 // Configuration will be linear.
3119 case 2:
3120 {
3121 // If we have no bonds yet of the 2 possible, we draw a site
3122 // along the y-axis and its inverse.
3123 if (bonds.size() == 0) {
3124 CPoint3D v(0.0f, 1.0f, 0.0f);
3125 DO_SITE(v, 1);
3126 if (lpCount < 1) {DO_SITE(v * -1.0f, 2);}
3127 }
3128
3129 else if (bonds.size() == 1) {
3130 DO_SITE(vecs[0] * -1.0f, 2);
3131 }
3132 }
3133 break;
3134
3135 // Configuration will be planar triangular.
3136 case 3:
3137 if (bonds.size() == 0) {
3138 CPoint3D v(0.0f, 1.0f, 0.0f);
3139 DO_SITE(v, 1);
3140 }
3141
3142 if (lpCount < 2) {
3143 CPoint3D v3;
3144 if (bonds.size() >= 2) { //Setup the 3rd vector to have an equal angle to both bonds
3145 //Just add the two bond vectors and invert
3146 // CPoint3D offset = lAtoms[secondaryBondedAtm].Position - lAtoms[iatom].Position;
3147 // Normalize3D(&offset);
3148 // v3 = b1Offset + offset;
3149 // Normalize3D(&v3);
3150 v3 = (vecs[0] + vecs[1]) * -1.0f;
3151 Normalize3D(&v3);
3152 } else { //Draw the 2nd position
3153 b = sqrt(3.0f / 4.0f);
3154 CPoint3D v2(-b, -0.5f, 0.0f), temp(b, -0.5f, 0.0f);
3155 Rotate3DOffset(rot,temp,&v3);
3156 temp = v2;
3157 Rotate3DOffset(rot,temp,&v2);
3158 //When dealing with a single bond check the bonded atom for bonds and align them as
3159 //appropriate
3160 if (bonds.size() == 1) {
3161 bool test = false;
3162 CPoint3D OtherBondVector;
3163 for (long i=0; i<NumBonds; i++) {
3164 if ((bonded_atoms[0] == lBonds[i].Atom1)&&(lBonds[i].Order > kHydrogenBond)&&
3165 (lBonds[i].Atom2 != iatom)) {
3166 OtherBondVector = lAtoms[lBonds[i].Atom2].Position - lAtoms[bonded_atoms[0]].Position;
3167 test = true;
3168 break;
3169 } else if ((bonded_atoms[0] == lBonds[i].Atom2)&&(lBonds[i].Order > kHydrogenBond)&&
3170 (lBonds[i].Atom1 != iatom)) {
3171 OtherBondVector = lAtoms[lBonds[i].Atom1].Position - lAtoms[bonded_atoms[0]].Position;
3172 test = true;
3173 break;
3174 }
3175 }
3176 if (test) {
3177 Normalize3D(&OtherBondVector);
3178 float dp = DotProduct3D(&OtherBondVector, &b1Offset);
3179 if (fabs(dp) < 0.99) {
3180 CPoint3D cross;
3181 CrossProduct3D(&b1Offset, &OtherBondVector, &cross);
3182 Normalize3D(&cross);
3183 CrossProduct3D(&b1Offset, &cross, &OtherBondVector);
3184 v2 = b1Offset*-0.5 + OtherBondVector*-0.866025;
3185 v3 = b1Offset*-0.5 + OtherBondVector*0.866025;
3186 }
3187 }
3188 }
3189 if (bonds.size() <= 1) {
3190 DO_SITE(v2, 2);
3191 }
3192 }
3193 if (lpCount < 1) {
3194 DO_SITE(v3, 3);
3195 }
3196 }
3197 break;
3198
3199 case 4:
3200 if (bonds.size() == 0) {
3201 CPoint3D v(0.0f, 1.0f, 0.0f);
3202 DO_SITE(v, 1);
3203 }
3204 if (lpCount < 3) {
3205
3206 // If three bonds are already formed, we calculate the fourth
3207 // as the sum of the three, inverted.
3208 if (bonds.size() == 3) {
3209 //obtain the bond vector by summing the bond vectors, normalize and invert
3210 CPoint3D sum(0.0,0.0,0.0), temp, bonds[3];
3211 int found=0;
3212 for (int i=0; i<NumBonds; i++) {
3213 if (lBonds[i].Atom1 == iatom) {
3214 temp = lAtoms[lBonds[i].Atom2].Position - lAtoms[iatom].Position;
3215 Normalize3D(&temp);
3216 bonds[found] = temp;
3217 found++;
3218 } else if (lBonds[i].Atom2 == iatom) {
3219 temp = lAtoms[lBonds[i].Atom1].Position - lAtoms[iatom].Position;
3220 Normalize3D(&temp);
3221 bonds[found] = temp;
3222 found++;
3223 }
3224 }
3225 CrossProduct3D(&(bonds[0]),&(bonds[1]),&sum);
3226 Normalize3D(&sum);
3227 float direction = DotProduct3D(&sum,&(bonds[2]));
3228 if (direction > 0.0) direction = -1.0;
3229 else direction = 1.0;
3230 CrossProduct3D(&(bonds[1]),&(bonds[2]),&temp);
3231 Normalize3D(&temp);
3232 sum += temp;
3233 CrossProduct3D(&(bonds[2]),&(bonds[0]),&temp);
3234 Normalize3D(&temp);
3235 sum += temp;
3236 Normalize3D(&sum);
3237 sum *= direction;
3238 DO_SITE(sum, 4);
3239 }
3240
3241 // If only one bond is already formed, we calculate the second, third
3242 // and fourth as ....
3243 else if (bonds.size() < 2) {
3244 c = -1.0f/3.0f; //cos(109.5)
3245 s = sqrt(1.0f-(c*c)); //sin(109.5)
3246 b = s * sin(kPi/3.0);
3247 CPoint3D v2, temp(0.0f, c, s), v3(b, c, -0.5*s), v4(-b, c, -0.5*s);
3248 Rotate3DOffset(rot,temp,&v2);
3249 DO_SITE(v2, 2);
3250 if (lpCount < 2) {
3251 temp = v3;
3252 Rotate3DOffset(rot,temp,&v3);
3253 DO_SITE(v3, 3);
3254 if (lpCount < 1) {
3255 temp = v4;
3256 Rotate3DOffset(rot,temp,&v4);
3257 DO_SITE(v4, 4);
3258 }
3259 }
3260 }
3261
3262 // If two bonds are already formed, we calculate the
3263 // remaining two as ...
3264 else {
3265 CPoint3D offset = lAtoms[bonded_atoms[1]].Position - lAtoms[iatom].Position;
3266 Normalize3D(&offset);
3267 CPoint3D cross;
3268 CrossProduct3D(&b1Offset,&offset,&cross);
3269 offset += b1Offset;
3270 offset *= -1.0f;
3271 Normalize3D(&offset);
3272 Normalize3D(&cross);
3273 CPoint3D v3t;
3274 float c = cos(0.5*acos(-1.0f/3.0f)); //cos(109.5/2)
3275 float s = sqrt(1.0f-(c*c)); //sin(109.5/2)
3276 v3t = offset*c + cross*s;
3277 Normalize3D(&v3t);
3278 DO_SITE(v3t, 3);
3279 if (lpCount == 0) {
3280 v3t = offset*c - cross*s;
3281 Normalize3D(&v3t);
3282 DO_SITE(v3t, 4);
3283 }
3284 }
3285 }
3286 break;
3287
3288 case 5:
3289 // Here axial and equitorial sites are different. LPs should go in
3290 // axial sites first.
3291 // Hmm... Unclear why the previous comment was made. LPs prefer to be as far from other
3292 // items as possible so should be in equitorial positions since axial are 90 degrees from
3293 // 3 sites while equitorial are 90 degrees from 2 sites and 120 degrees from 2 sites.
3294
3295 if (bonds.size() == 0) {
3296 float c = sqrtf(3) / 2.0f;
3297 CPoint3D axial(0.0f, 1.0f, 0.0f);
3298 CPoint3D equat1(0.0f, 0.0f, -1.0f);
3299 CPoint3D equat2(c, 0.0f, 0.5f);
3300 CPoint3D equat3(-c, 0.0f, 0.5f);
3301
3302 DO_SITE(axial, 4);
3303 if (lpCount < 4) {
3304 DO_SITE(axial * -1.0f, 5);
3305 if (lpCount < 3) {
3306 DO_SITE(equat1, 1);
3307 if (lpCount < 2) {
3308 DO_SITE(equat2, 2);
3309 if (lpCount < 1) {DO_SITE(equat3, 3);}
3310 }
3311 }
3312 }
3313 }
3314
3315 else if (bonds.size() == 1) {
3316
3317 CPoint3D equat1;
3318 CPoint3D equat2;
3319 CPoint3D equat3;
3320 Matrix4D rotate_mat;
3321
3322 OrthoVector(vecs[0], equat1);
3323 RotateAroundAxis(rotate_mat, vecs[0], 120.0f);
3324 Rotate3DPt(rotate_mat, equat1, &equat2);
3325 Normalize3D(&equat2);
3326 equat3 = (equat1 + equat2) * -1.0f;
3327 Normalize3D(&equat3);
3328
3329 DO_SITE(vecs[0] * -1.0f, 5);
3330 if (lpCount < 3) {
3331 DO_SITE(equat1, 1);
3332 if (lpCount < 2) {
3333 DO_SITE(equat2, 2);
3334 if (lpCount < 1) {DO_SITE(equat3, 3);}
3335 }
3336 }
3337 }
3338
3339 else if (bonds.size() == 2) {
3340
3341 float dot = DotProduct3D(&vecs[0], &vecs[1]);
3342
3343 // If the two vectors are close to orthogonal, we assume the
3344 // primary bond is axial and the secondary equitorial.
3345 if (fabs(dot) < 0.25f) {
3346 Matrix4D rotate_mat;
3347 CPoint3D equat2;
3348 CPoint3D equat3;
3349
3350 if (lpCount < 3) {
3351 DO_SITE(vecs[0] * -1.0f, 5);
3352 if (lpCount < 2) {
3353 RotateAroundAxis(rotate_mat, vecs[0], 120.0f);
3354 Rotate3DPt(rotate_mat, vecs[1], &equat2);
3355 Normalize3D(&equat2);
3356 equat3 = (vecs[1] + equat2) * -1.0f;
3357 Normalize3D(&equat3);
3358
3359 DO_SITE(equat2, 2);
3360 if (lpCount < 1) {DO_SITE(equat3, 3);}
3361 }
3362 }
3363 } else {
3364 if (lpCount < 3) {
3365 CPoint3D equat1;
3366 CPoint3D equat2;
3367 CPoint3D equat3;
3368 Matrix4D rotate_mat;
3369
3370 OrthoVector(vecs[0], equat1);
3371 RotateAroundAxis(rotate_mat, vecs[0], 120.0f);
3372 Rotate3DPt(rotate_mat, equat1, &equat2);
3373 Normalize3D(&equat2);
3374 equat3 = (equat1 + equat2) * -1.0f;
3375 Normalize3D(&equat3);
3376
3377 DO_SITE(equat1, 1);
3378 if (lpCount < 2) {
3379 DO_SITE(equat2, 2);
3380 if (lpCount < 1) {DO_SITE(equat3, 3);}
3381 }
3382 }
3383 }
3384 }
3385
3386 else if (bonds.size() == 3) {
3387
3388 float dot12 = DotProduct3D(&vecs[0], &vecs[1]);
3389 float dot13 = DotProduct3D(&vecs[0], &vecs[2]);
3390 float dot23 = DotProduct3D(&vecs[1], &vecs[2]);
3391
3392 // Two vectors are opposite, meaning two are axial.
3393 if (dot12 < -0.75 || dot13 < -0.75f || dot23 < -0.75f) {
3394
3395 CPoint3D cross1;
3396 CPoint3D cross2;
3397 CPoint3D equat2;
3398 CPoint3D equat3;
3399 int axial_id=0;
3400 int equat_id=1;
3401 Matrix4D rotate_mat;
3402
3403 if (dot12 < -0.75) {
3404 axial_id = 0;
3405 equat_id = 2;
3406 } else if (dot13 < -0.75f) {
3407 axial_id = 0;
3408 equat_id = 1;
3409 } else if (dot23 < -0.75f) {
3410 axial_id = 1;
3411 equat_id = 0;
3412 }
3413
3414 CrossProduct3D(&vecs[axial_id], &vecs[equat_id], &cross1);
3415 CrossProduct3D(&cross1, &vecs[equat_id], &cross2);
3416 Normalize3D(&cross2);
3417
3418 RotateAroundAxis(rotate_mat, cross2, 120.0f);
3419 Rotate3DPt(rotate_mat, vecs[equat_id], &equat2);
3420 equat3 = (equat2 + vecs[equat_id]) * -1.0f;
3421 // Don't need to normalize equat3 since it is by
3422 // construction.
3423
3424 DO_SITE(equat2, 4);
3425 if (lpCount < 1) {DO_SITE(equat3, 5);}
3426 }
3427
3428 else if (fabs(dot12) < 0.25f || fabs(dot13) < 0.25f ||
3429 fabs(dot23) < 0.25f) {
3430
3431 int axial_id=0;
3432 int equat1_id=1;
3433 int equat2_id=2;
3434
3435 if (fabs(dot12) < 0.25) {
3436 equat1_id = 2;
3437 if (fabs(dot13) < 0.25f) {
3438 axial_id = 0;
3439 equat2_id = 1;
3440 } else {
3441 axial_id = 1;
3442 equat2_id = 0;
3443 }
3444 } else if (fabs(dot13) < 0.25f) {
3445 equat1_id = 1;
3446 if (fabs(dot12) < 0.25f) {
3447 axial_id = 0;
3448 equat2_id = 2;
3449 } else {
3450 axial_id = 2;
3451 equat2_id = 0;
3452 }
3453 } else if (fabs(dot23) < 0.25f) {
3454 equat1_id = 0;
3455 if (fabs(dot13) < 0.25f) {
3456 axial_id = 2;
3457 equat2_id = 1;
3458 } else {
3459 axial_id = 1;
3460 equat2_id = 2;
3461 }
3462 }
3463
3464 CPoint3D cross;
3465 CPoint3D equat3;
3466
3467 CrossProduct3D(&vecs[equat1_id], &vecs[equat2_id], &cross);
3468 Normalize3D(&cross);
3469 if (DotProduct3D(&cross, &vecs[axial_id]) > 0.0f) {
3470 cross *= -1.0f;
3471 }
3472
3473 equat3 = (vecs[equat1_id] + vecs[equat2_id]) * -1.0f;
3474 Normalize3D(&equat3);
3475
3476 DO_SITE(cross, 5);
3477 if (lpCount < 1) {DO_SITE(equat3, 4);}
3478 }
3479
3480 else {
3481 CPoint3D cross;
3482
3483 CrossProduct3D(&vecs[0], &vecs[1], &cross);
3484 Normalize3D(&cross);
3485
3486 DO_SITE(cross, 4);
3487 if (lpCount < 1) {DO_SITE(cross * -1.0f, 5);}
3488 }
3489
3490 }
3491
3492 else if (bonds.size() == 4) {
3493
3494 float dot_ij;
3495 bool is_axial = false;
3496 bool is_equatorial;
3497 int i, j;
3498
3499 for (i = 0; i < 4; i++) {
3500 is_equatorial = false;
3501 for (j = 0; j < 4; j++) {
3502 if (i == j) continue;
3503
3504 dot_ij = DotProduct3D(&vecs[i], &vecs[j]);
3505
3506 // Vector has an antiparallel mate, so we can't stop
3507 // now.
3508 if (dot_ij < -0.75) {
3509 is_axial = true;
3510 break;
3511 }
3512
3513 // Vector has a non-orthogonal mate (which must also
3514 // be non-antiparallel if we got here), so it can't
3515 // be an axial vector.
3516 else if (fabs(dot_ij) > 0.25) {
3517 is_equatorial = true;
3518 break;
3519 }
3520 }
3521
3522 if (is_axial) {
3523 CPoint3D equat3;
3524 int ids[4] = {0, 1, 2, 3};
3525 int tmp;
3526 tmp = ids[0]; ids[0] = i; ids[i] = tmp;
3527 tmp = ids[1]; ids[1] = j; ids[j] = tmp;
3528 equat3 = (vecs[ids[2]] + vecs[ids[3]]) * -1.0f;
3529 Normalize3D(&equat3);
3530 DO_SITE(equat3, 5);
3531 break;
3532 } else if (!is_equatorial) {
3533 DO_SITE(vecs[i] * -1.0f, 5);
3534 break;
3535 }
3536 }
3537
3538 if (!is_axial && is_equatorial) {
3539 CPoint3D cross;
3540 CrossProduct3D(&vecs[0], &vecs[1], &cross);
3541 Normalize3D(&cross);
3542 DO_SITE(cross, 5);
3543 }
3544 }
3545
3546 break;
3547 case 6:
3548
3549 if (bonds.size() == 0) {
3550 CPoint3D e1(1.0f, 0.0f, 0.0f);
3551 CPoint3D e2(0.0f, 1.0f, 0.0f);
3552 CPoint3D e3(0.0f, 0.0f, 1.0f);
3553
3554 DO_SITE(e1, 1);
3555 if (lpCount < 5) {
3556 DO_SITE(e1 * -1.0f, 2);
3557 if (lpCount < 4) {
3558 DO_SITE(e2, 3);
3559 if (lpCount < 3) {
3560 DO_SITE(e2 * -1.0f, 4);
3561 if (lpCount < 2) {
3562 DO_SITE(e3, 5);
3563 if (lpCount < 1) {DO_SITE(e3 * -1.0f, 6);}
3564 }
3565 }
3566 }
3567 }
3568 }
3569
3570 else if (bonds.size() == 1) {
3571 CPoint3D ortho;
3572 CPoint3D cross;
3573
3574 OrthoVector(vecs[0], ortho);
3575 CrossProduct3D(&vecs[0], &ortho, &cross);
3576
3577 DO_SITE(vecs[0] * -1.0f, 2);
3578 if (lpCount < 4) {
3579 DO_SITE(ortho, 3);
3580 if (lpCount < 3) {
3581 DO_SITE(ortho * -1.0f, 4);
3582 if (lpCount < 2) {
3583 DO_SITE(cross, 5);
3584 if (lpCount < 1) {DO_SITE(cross * -1.0f, 6);}
3585 }
3586 }
3587 }
3588 }
3589
3590 else if (bonds.size() == 2) {
3591
3592 // Use their dot product to determine if they're orthogonal or
3593 // otherwise.
3594 float dot = DotProduct3D(&vecs[0], &vecs[1]);
3595
3596 // If the two existing bonds are closer to being orthogonal
3597 // than parallel (or antiparallel), then we have four bonds
3598 // to compute. One bond will be the inverse of the primary
3599 // bond, two will be the cross product and its inverse of the
3600 // two bonds, and the fourth will be another cross product,
3601 // this time using the primary bond and the first cross
3602 // product. It might make more sense to use the inverse of
3603 // the secondary bond instead, but doing it this way makes at
3604 // least 5 of the 6 bonds octahedrally aligned.
3605 if (fabs(dot) < 0.7f) {
3606 CPoint3D cross;
3607 CPoint3D cross2;
3608 CrossProduct3D(&vecs[0], &vecs[1], &cross);
3609 Normalize3D(&cross);
3610 CrossProduct3D(&vecs[0], &cross, &cross2);
3611 Normalize3D(&cross2);
3612
3613 if (DotProduct3D(&cross2, &vecs[1]) > 0.0f) {
3614 cross2 *= -1.0f;
3615 }
3616
3617 DO_SITE(vecs[0] * -1.0f, 3);
3618 if (lpCount < 3) {
3619 DO_SITE(cross2, 4);
3620 if (lpCount < 2) {
3621 DO_SITE(cross, 5);
3622 if (lpCount < 1) {DO_SITE(cross * -1.0f, 6);}
3623 }
3624 }
3625 }
3626
3627 // Otherwise, the two existing bonds must be antiparallel.
3628 // In that case, we calculate the other four bonds like so:
3629 // one bond is chosen from the many that are orthogonal to
3630 // the primary bond, the next is its inverse, and the last
3631 // two are the crossproduct and its inverse of the primary
3632 // bond and its orthogonal mate.
3633 else {
3634 CPoint3D ortho;
3635 CPoint3D cross;
3636
3637 OrthoVector(vecs[0], ortho);
3638 CrossProduct3D(&vecs[0], &ortho, &cross);
3639 Normalize3D(&cross);
3640
3641 DO_SITE(ortho, 3);
3642 if (lpCount < 3) {
3643 DO_SITE(ortho * -1.0f, 4);
3644 if (lpCount < 2) {
3645 DO_SITE(cross, 5);
3646 if (lpCount < 1) {DO_SITE(cross * -1.0f, 6);}
3647 }
3648 }
3649 }
3650
3651 }
3652
3653 else if (bonds.size() == 3) {
3654
3655 float dot12 = DotProduct3D(&vecs[0], &vecs[1]);
3656 float dot13 = DotProduct3D(&vecs[0], &vecs[2]);
3657 float dot23 = DotProduct3D(&vecs[1], &vecs[2]);
3658
3659 // If all three vectors are close to pairwise orthogonal, then
3660 // they form a basis. Unfortunately, they may not be exactly
3661 // independent. Since we want to suggest orthogonal sites, we
3662 // choose the primary vector's inverse as the first site. The
3663 // second is calculated as the cross product of the primary
3664 // and secondary bonds, and the third is the cross product of
3665 // the primary bond and the first cross product. This makes
3666 // all three orthogonal.
3667 if (fabs(dot12) < 0.7f && fabs(dot13) < 0.7f &&
3668 fabs(dot23) < 0.7f) {
3669
3670 CPoint3D cross1, cross2;
3671 CrossProduct3D(&vecs[0], &vecs[1], &cross1);
3672 Normalize3D(&cross1);
3673 CrossProduct3D(&vecs[0], &cross1, &cross2);
3674 Normalize3D(&cross2);
3675
3676 // The cross products should point in the opposite
3677 // direction of the existing bonds, so we may need to
3678 // flip them since they can end up coincident with them.
3679 if (DotProduct3D(&cross1, &vecs[2]) > 0.0f) {
3680 cross1 = cross1 * -1.0f;
3681 }
3682
3683 if (DotProduct3D(&cross2, &vecs[1]) > 0.0f) {
3684 cross2 = cross2 * -1.0f;
3685 }
3686
3687 DO_SITE(vecs[0] * -1.0f, 4);
3688 if (lpCount < 2) {
3689 DO_SITE(cross1, 5);
3690 if (lpCount < 1) {DO_SITE(cross2, 6);}
3691 }
3692 }
3693
3694 // Otherwise, two of the three bonds are antiparallel. In this
3695 // case, we figure out which bond doesn't have an inverse and
3696 // one that does. The first suggested bond is the inverse of
3697 // the unopposed bond, and the second and third are the cross
3698 // product and its inverse of the unpaired bond and one of the
3699 // paired bonds.
3700 else {
3701 CPoint3D paired_vec;
3702 CPoint3D ortho_vec;
3703
3704 if (dot12 < -0.7f) {
3705 paired_vec = vecs[0];
3706 ortho_vec = vecs[2];
3707 } else if (dot13 < -0.7f) {
3708 paired_vec = vecs[0];
3709 ortho_vec = vecs[1];
3710 } else {
3711 paired_vec = vecs[2];
3712 ortho_vec = vecs[0];
3713 }
3714
3715 CPoint3D cross;
3716 CrossProduct3D(&paired_vec, &ortho_vec, &cross);
3717 Normalize3D(&cross);
3718
3719 DO_SITE(ortho_vec * -1.0f, 4);
3720 if (lpCount < 2) {
3721 DO_SITE(cross, 5);
3722 if (lpCount < 1) {DO_SITE(cross * -1.0f, 6);}
3723 }
3724 }
3725
3726 }
3727
3728 else if (bonds.size() == 4) {
3729
3730 int i, j;
3731 int nuninverted = 0;
3732 int uninverted[2];
3733 int ortho1 = -1, ortho2 = -1;
3734 bool has_no_inverse;
3735
3736 for (i = 0; i < 4; i++) {
3737 has_no_inverse = true;
3738 for (j = 0; j < 4; j++) {
3739 if (i != j) {
3740 if (DotProduct3D(&vecs[i], &vecs[j]) < -0.7f) {
3741 has_no_inverse = false;
3742 break;
3743 } else {
3744 ortho1 = i;
3745 ortho2 = j;
3746 }
3747 }
3748 }
3749 if (has_no_inverse) {
3750 uninverted[nuninverted] = i;
3751 nuninverted++;
3752 }
3753 }
3754
3755 if (nuninverted == 2) {
3756 DO_SITE(vecs[uninverted[0]] * -1.0f, 5);
3757 if (lpCount == 0) {DO_SITE(vecs[uninverted[1]] * -1.0f, 6);}
3758 }
3759
3760 else {
3761 CPoint3D cross;
3762 CrossProduct3D(&vecs[ortho1], &vecs[ortho2], &cross);
3763 Normalize3D(&cross);
3764
3765 DO_SITE(cross, 5);
3766 if (lpCount == 0) {DO_SITE(cross * -1.0f, 6);}
3767 }
3768
3769 }
3770
3771 else if (bonds.size() == 5) {
3772
3773 int i, j;
3774 bool has_no_inverse;
3775
3776 // If there isn't an uninverted bond, just go opposite the
3777 // primary bond.
3778 int uninverted = 0;
3779
3780 for (i = 0; i < 5; i++) {
3781 has_no_inverse = true;
3782 for (j = 0; j < 5; j++) {
3783 if (i != j) {
3784 if (DotProduct3D(&vecs[i], &vecs[j]) < -0.7f) {
3785 has_no_inverse = false;
3786 break;
3787 }
3788 }
3789 }
3790 if (has_no_inverse) {
3791 uninverted = i;
3792 break;
3793 }
3794 }
3795
3796 DO_SITE(vecs[uninverted] * -1.0f, 6);
3797 }
3798
3799 break;
3800 }
3801
3802 if (site_id == 0) {
3803 glPopName();
3804 }
3805
3806 }
3807
3808 /* ------------------------------------------------------------------------- */
3809
DrawBond(const Bond & bond,const mpAtom & atom1,const mpAtom & atom2,const WinPrefs & Prefs,GLUquadric * quadric,GLdouble * modelview,GLdouble * proj,GLint * viewport,GLuint sphere_list,bool highlighting_on,bool cap_dependent)3810 void MolDisplayWin::DrawBond(const Bond& bond, const mpAtom& atom1,
3811 const mpAtom& atom2, const WinPrefs& Prefs,
3812 GLUquadric *quadric, GLdouble *modelview,
3813 GLdouble *proj, GLint *viewport,
3814 GLuint sphere_list, bool highlighting_on,
3815 bool cap_dependent) {
3816
3817 CPoint3D v1, v2, offset, NormalOffset, NormEnd;
3818 CPoint3D NormStart(0, 0, 1);
3819 Matrix4D rotMat;
3820 long Quality = (long) Prefs.GetQD3DAtomQuality();
3821 float AtomScale = Prefs.GetAtomScale();
3822 GLdouble BondSize = Prefs.GetQD3DBondWidth() * 0.5;
3823
3824 BondOrder logical_order = bond.Order;
3825
3826 if (logical_order == kAromaticBond) logical_order = kDoubleBond;
3827 if (!Prefs.ColorBondHalves()) logical_order = kSingleBond; //only generate multiple pipes when colored by atom color
3828
3829 GLdouble tmpBondSize = BondSize / MAX(logical_order, 1);
3830 if (logical_order > 1) tmpBondSize *= 1.5;
3831 GLdouble baseBondOffset = -1.75 * tmpBondSize * (MAX(logical_order, 1) - 1);
3832
3833 if (atom1.GetInvisibility() || atom2.GetInvisibility()) return;
3834
3835 // If this is a bond that's going to need multiple pipes, we want each pipe
3836 // to be shifted in screen space so that all pipes are visible all the
3837 // time.
3838 CPoint3D offset_vec; // Direction to shift bond cylinders
3839 if (logical_order > 1) {
3840 GLdouble scr_coords1[3]; // Screen coordinates of atom1
3841 GLdouble scr_coords2[3]; // Screen coordinates of atom2
3842 CPoint3D scr_vec; // Screen space vector between atoms
3843 GLdouble perp_obj[3]; // Object coords on vector perp. to scr_vec
3844
3845 // Find screen coordinates of one atom.
3846 gluProject(atom1.Position.x, atom1.Position.y, atom1.Position.z,
3847 modelview, proj, viewport,
3848 &(scr_coords1[0]), &(scr_coords1[1]), &(scr_coords1[2]));
3849
3850 // Find screen coordinates of other atom.
3851 gluProject(atom2.Position.x, atom2.Position.y, atom2.Position.z,
3852 modelview, proj, viewport,
3853 &(scr_coords2[0]), &(scr_coords2[1]), &(scr_coords2[2]));
3854
3855 // Find vector perpendicular to vector between two screen points and
3856 // normalize it so we can scalar multiply it later. We flip and
3857 // negate the slope of the line between the two screen coordinates to
3858 // get the slop of the perpendicular line.
3859 scr_vec.x = scr_coords2[1] - scr_coords1[1];
3860 scr_vec.y = scr_coords1[0] - scr_coords2[0];
3861 scr_vec.z = 0;
3862 scr_vec *= 1.0f / scr_vec.Magnitude();
3863
3864 // Now find a point on the perpendicular vector with atom1's depth
3865 // and get its object coordinates.
3866 gluUnProject(scr_coords1[0] + scr_vec.x * 10.0f,
3867 scr_coords1[1] + scr_vec.y * 10.0f,
3868 scr_coords1[2],
3869 modelview, proj, viewport,
3870 &(perp_obj[0]), &(perp_obj[1]), &(perp_obj[2]));
3871
3872 // Finally, we see what direction all bond cylinders must be offset
3873 // so that they will always stay in view.
3874 offset_vec = CPoint3D(perp_obj[0], perp_obj[1], perp_obj[2]) -
3875 atom1.Position;
3876 offset_vec *= 1.0f / offset_vec.Magnitude();
3877 }
3878
3879 // If there's only one pipe, we need no offset.
3880 else {
3881 offset_vec = CPoint3D(0.0f, 0.0f, 0.0f);
3882 }
3883
3884 // For each "sub-bond" between these two atoms...
3885 for (int ipipe = 0; ipipe < MAX(logical_order, 1); ++ipipe) {
3886
3887 v1 = atom1.Position + offset_vec * baseBondOffset +
3888 offset_vec * tmpBondSize * ipipe * 3.5;
3889 v2 = atom2.Position + offset_vec * baseBondOffset +
3890 offset_vec * tmpBondSize * ipipe * 3.5;
3891
3892 offset = v2 - v1;
3893 float length = offset.Magnitude();
3894 if (length > 0.00001) {
3895 NormalOffset = offset * (1.0f / length);
3896 } else {
3897 NormalOffset = CPoint3D(0, 0, 0);
3898 }
3899
3900 SetRotationMatrix(rotMat, &NormStart, &NormalOffset);
3901
3902 rotMat[3][0] = v1.x;
3903 rotMat[3][1] = v1.y;
3904 rotMat[3][2] = v1.z;
3905
3906 glPushMatrix(); // bond pipe
3907 glMultMatrixf((const GLfloat *) &rotMat);
3908
3909 // Both hydrogen bonds and one part of the aromatic bond and drawn
3910 // as a line of spheres.
3911 if (bond.Order == kHydrogenBond ||
3912 (bond.Order == kAromaticBond && ipipe == 1)) {
3913 DrawPipeSpheres(Prefs, length, BondSize,
3914 highlighting_on && !bond.GetSelectState(),
3915 sphere_list);
3916 glPopMatrix(); // bond pipe
3917 continue;
3918 }
3919
3920 // Now, if a bond is selected but not this one, we need to draw an
3921 // encapsulating cylinder to mask it out.
3922 if (highlighting_on && !bond.GetSelectState()) {
3923
3924 glEnable(GL_POLYGON_STIPPLE);
3925 glPolygonStipple(stippleMask);
3926
3927 glColor3f(0.0f, 0.0f, 0.0f);
3928 DrawPipeCylinder(length, quadric, Prefs.DrawWireFrame() ? 3 : 0,
3929 sphere_list, tmpBondSize * 1.01f, Quality);
3930
3931 glDisable(GL_POLYGON_STIPPLE);
3932
3933 }
3934
3935 // We may need to draw two cylinders if the user wants the bonds
3936 // colored according to their element.
3937 if (Prefs.ColorBondHalves()) {
3938 // Center the color change at the middle of the visible part of the
3939 // bond.
3940 float radius1 = AtomScale*Prefs.GetAtomSize(atom1.GetType() - 1);
3941 float radius2 = AtomScale*Prefs.GetAtomSize(atom2.GetType() - 1);
3942 float percent1 = radius1/length;
3943 float percent2 = radius2/length;
3944 float centerPercent = 0.5 + 0.5 * (percent1 - percent2);
3945
3946 // Draw first half.
3947 Prefs.ChangeColorAtomColor(atom1.GetType());
3948 DrawPipeCylinder(length * centerPercent, quadric,
3949 (Prefs.DrawWireFrame() ||
3950 (cap_dependent && !atom1.IsSymmetryUnique()) ||
3951 (Prefs.ShowEFPWireFrame() &&
3952 atom1.IsEffectiveFragment())) ? 1 : 0,
3953 sphere_list,
3954 tmpBondSize, Quality);
3955
3956 // Draw second half starting at the weighted halfway point.
3957 glPushMatrix(); // halfway point
3958 glTranslatef(0.0f, 0.0f, length * centerPercent);
3959 Prefs.ChangeColorAtomColor(atom2.GetType());
3960 DrawPipeCylinder(length * (1.0f - centerPercent), quadric,
3961 (Prefs.DrawWireFrame() ||
3962 (cap_dependent && !atom2.IsSymmetryUnique()) ||
3963 (Prefs.ShowEFPWireFrame() &&
3964 atom2.IsEffectiveFragment())) ? 2 : 0,
3965 sphere_list, tmpBondSize, Quality);
3966 glPopMatrix(); // halfway point
3967
3968 }
3969
3970 // We only need to draw one cylinder the whole length of a bond since
3971 // the user's not interested in coloring each half differently.
3972 else {
3973 Prefs.ChangeColorBondColor(bond.Order);
3974 DrawPipeCylinder(length, quadric,
3975 (Prefs.DrawWireFrame() ||
3976 (cap_dependent &&
3977 (!atom1.IsSymmetryUnique() ||
3978 !atom2.IsSymmetryUnique())) ||
3979 (Prefs.ShowEFPWireFrame() &&
3980 (atom1.IsEffectiveFragment() ||
3981 atom2.IsEffectiveFragment()))) ? 3 : 0,
3982 sphere_list, tmpBondSize, Quality);
3983 }
3984
3985 glPopMatrix(); // bond pipe
3986 }
3987
3988 }
3989
3990 /* ------------------------------------------------------------------------- */
3991 /**
3992 * This function draws a series of spheres from the origin along the negative
3993 * z-axis for length units. The spheres are optionally drawn with a hatching
3994 * mask.
3995 * @param Prefs Prefences used to determine sphere color.
3996 * @param length Distance along z-axis along which to draw spheres.
3997 * @param scale_factor Size of spheres relative to current transform.
3998 * @param is_masked Flag indicating whether to draw stippled overlay to mask
3999 * spheres.
4000 * @param sphere_list Preassembled display list that draws a sphere.
4001 */
4002
DrawPipeSpheres(const WinPrefs & Prefs,float length,float scale_factor,bool is_masked,GLuint sphere_list)4003 void DrawPipeSpheres(const WinPrefs& Prefs, float length, float scale_factor,
4004 bool is_masked, GLuint sphere_list) {
4005
4006 glPushMatrix(); // offset from origin
4007
4008 // Plot as a series of spheres, starting slightly offset from the origin.
4009 float pos = 0.75f * scale_factor;
4010 glTranslatef(0.0, 0.0, pos);
4011 while (pos < length) {
4012 glPushMatrix(); // scaled sphere location
4013 glScalef(scale_factor, scale_factor, scale_factor);
4014
4015 Prefs.ChangeColorBondColor(kHydrogenBond);
4016 glCallList(sphere_list);
4017
4018 if (is_masked) {
4019 glEnable(GL_POLYGON_STIPPLE);
4020 glPolygonStipple(stippleMask);
4021 glPushMatrix(); // stipple-scaled sphere location
4022 glScalef(1.01f, 1.01f, 1.01f);
4023 glColor3f(0.0f, 0.0f, 0.0f);
4024 glCallList(sphere_list);
4025 glPopMatrix(); // stipple-scaled sphere location
4026 glDisable(GL_POLYGON_STIPPLE);
4027 }
4028
4029 glPopMatrix(); // scaled sphere location
4030
4031 glTranslatef(0.0f, 0.0f, 2.5f * scale_factor);
4032 pos += 2.5f * scale_factor;
4033 }
4034
4035 glPopMatrix(); // offset from origin
4036
4037 }
4038
4039 /* ------------------------------------------------------------------------- */
4040 /**
4041 * Draws a cylinder from the current transformation's origin along negative z
4042 * axis for length units. Assumes modeling transformation is in place, and
4043 * OpenGL commands are ready to be issued. Useful for drawing bond pipes that
4044 * are optionally capped.
4045 * @param length Extent of the cylinder along the z-axis.
4046 * @param quadric Preallocated GLU quadric object.
4047 * @param ncaps Number of end caps to draw. Should be 1 to draw endcap only at
4048 * origin, 2 only at the opposite end, and 3 to draw both.
4049 * @param sphere_list Preassembled display list id for drawing endcaps, which
4050 * should be valid if ncaps != 0.
4051 * @param radius Radius of cylinder.
4052 * @param quality Measure of "facetedness" of cylinder's geometric
4053 * representation.
4054 */
4055
DrawPipeCylinder(float length,GLUquadric * quadric,unsigned int ncaps,GLuint sphere_list,float radius,long quality)4056 void DrawPipeCylinder(float length, GLUquadric *quadric, unsigned int ncaps,
4057 GLuint sphere_list, float radius, long quality) {
4058
4059 gluCylinder(quadric, radius, radius, length, quality,
4060 (long) (0.5f * quality));
4061
4062 // Draw sphere at end0 if requested.
4063 if (ncaps & 1) {
4064 glPushMatrix(); // end0
4065 glScalef(radius, radius, radius);
4066 glCallList(sphere_list);
4067 glPopMatrix(); // end0
4068 }
4069
4070 // Draw sphere at end1 if requested.
4071 if (ncaps & 2) {
4072 glPushMatrix(); // end1
4073 glTranslatef(0.0f, 0.0f, length);
4074 glScalef(radius, radius, radius);
4075 glCallList(sphere_list);
4076 glPopMatrix(); // end1
4077 }
4078
4079 }
4080
4081 /* ------------------------------------------------------------------------- */
4082
GetShaderProgramFromFiles(const std::string & vert_filename,const std::string & frag_filename)4083 GLuint GetShaderProgramFromFiles(const std::string& vert_filename,
4084 const std::string& frag_filename) {
4085
4086 std::string vert_src;
4087 std::string frag_src;
4088
4089 if (!vert_filename.empty()) {
4090 if (!FileToString(vert_filename, vert_src)) {
4091 return 0;
4092 }
4093 }
4094
4095 if (!frag_filename.empty()) {
4096 if (!FileToString(frag_filename, frag_src)) {
4097 return 0;
4098 }
4099 }
4100
4101 return GetShaderProgram(vert_src, frag_src);
4102
4103 }
4104
4105 /* ------------------------------------------------------------------------- */
4106
GetShaderProgram(const std::string & vert_src,const std::string & frag_src)4107 GLuint GetShaderProgram(const std::string& vert_src,
4108 const std::string& frag_src) {
4109
4110 GLuint shader_prog = 0;
4111
4112 GLint succeeded;
4113 char *log;
4114 GLint log_length=0;
4115 GLint chars_written;
4116
4117 shader_prog = glCreateProgram();
4118
4119 if (vert_src.length()) {
4120 GLuint vert_shader;
4121 const char *v2 = &vert_src.c_str()[0];
4122
4123 // Vertex shader.
4124 vert_shader = glCreateShader(GL_VERTEX_SHADER);
4125 glShaderSource(vert_shader, 1, &v2, NULL);
4126 glCompileShader(vert_shader);
4127 glGetShaderiv(vert_shader, GL_COMPILE_STATUS, &succeeded);
4128
4129 glGetShaderiv(vert_shader, GL_INFO_LOG_LENGTH, &log_length);
4130 if (log_length > 1) {
4131 log = new char[log_length];
4132 glGetShaderInfoLog(vert_shader, log_length, &chars_written, log);
4133 std::cout << "log: " << log << std::endl;
4134 delete[] log;
4135 }
4136
4137 if (succeeded != GL_TRUE) {
4138 std::cout << "Vertex Program: " << vert_src << std::endl;
4139 wxLogMessage(wxT("Something went wrong with the shader."));
4140 }
4141
4142 glAttachShader(shader_prog, vert_shader);
4143 }
4144
4145 // Fragment shader.
4146 if (frag_src.length()) {
4147 GLuint frag_shader;
4148 const char *f2 = &frag_src.c_str()[0];
4149 frag_shader = glCreateShader(GL_FRAGMENT_SHADER);
4150 glShaderSource(frag_shader, 1, &f2, NULL);
4151
4152 glCompileShader(frag_shader);
4153 glGetShaderiv(frag_shader, GL_COMPILE_STATUS, &succeeded);
4154
4155 glGetShaderiv(frag_shader, GL_INFO_LOG_LENGTH, &log_length);
4156 if (log_length > 1) {
4157 log = new char[log_length];
4158 glGetShaderInfoLog(frag_shader, log_length, &chars_written, log);
4159 std::cout << "log: " << log << std::endl;
4160 delete[] log;
4161 }
4162
4163 if (succeeded != GL_TRUE) {
4164 std::cout << "Fragment Program: " << frag_src << std::endl;
4165 wxLogMessage(wxT("Something went wrong with the shader."));
4166 }
4167
4168 glAttachShader(shader_prog, frag_shader);
4169 }
4170
4171 glLinkProgram(shader_prog);
4172 glGetProgramiv(shader_prog, GL_LINK_STATUS, &succeeded);
4173 if (succeeded != GL_TRUE) {
4174 char *log;
4175 GLint log_length;
4176 GLint chars_written;
4177 glGetProgramiv(shader_prog, GL_INFO_LOG_LENGTH, &log_length);
4178 if (log_length) {
4179 log = new char[log_length];
4180 glGetProgramInfoLog(shader_prog, log_length, &chars_written, log);
4181 std::cout << "log: " << log << std::endl;
4182 delete[] log;
4183 exit(1);
4184 } else {
4185 wxLogMessage(wxT("Something went wrong with the shader program."));
4186 }
4187 }
4188
4189 return shader_prog;
4190
4191 }
4192
4193 /* ------------------------------------------------------------------------- */
4194
4195 /**
4196 This function (re)generates any structures that are dependent on preference
4197 settings. It should be called when preferences change. It assumes that a valid
4198 OpenGL context exists and is current, so it should not be called by a parent's
4199 constructor.
4200 */
DoPrefDependent()4201 void MolDisplayWin::DoPrefDependent() {
4202
4203 // All atoms are drawn using one display list that draws a sphere.
4204 // We delete any current list and create it anew.
4205 if (sphere_list) {
4206 glDeleteLists(sphere_list, 1);
4207 }
4208
4209 GLUquadric *quad = gluNewQuadric();
4210 gluQuadricOrientation(quad, GLU_OUTSIDE);
4211 gluQuadricNormals(quad, GLU_SMOOTH);
4212
4213 sphere_list = glGenLists(1);
4214 glNewList(sphere_list, GL_COMPILE);
4215
4216 #if 1
4217 float *verts;
4218 int *faces;
4219 int nverts;
4220 int nfaces;
4221 float *normals;
4222 int nlevels;
4223
4224 // We do some funny mapping from display quality (2 - 40) to number of
4225 // subdivisions. 2 maps to 0 levels, while all other values are mapped
4226 // such that 13 -> 3 levels and 40 -> 7 levels. Levels >7 recurse too
4227 // much.
4228 if (Prefs->GetQD3DAtomQuality() == 2) {
4229 nlevels = 0;
4230 } else {
4231 nlevels = (int) ((Prefs->GetQD3DAtomQuality() - 13.0f) *
4232 (4.0f / 27.0f) + 3.0f);
4233 }
4234
4235 GenerateOctahedron(nlevels, &verts, nverts, &faces, nfaces, &normals);
4236 glVertexPointer(3, GL_FLOAT, 0, verts);
4237
4238 glEnableClientState(GL_VERTEX_ARRAY);
4239 glEnableClientState(GL_NORMAL_ARRAY);
4240 glVertexPointer(3, GL_FLOAT, 0, verts);
4241 glNormalPointer(GL_FLOAT, 0, normals);
4242 glDrawElements(GL_TRIANGLES, nfaces * 3, GL_UNSIGNED_INT, faces);
4243 glDisableClientState(GL_NORMAL_ARRAY);
4244 glDisableClientState(GL_VERTEX_ARRAY);
4245
4246 delete[] verts;
4247 delete[] faces;
4248 delete[] normals;
4249 #else
4250 gluSphere(quad, 1.0f, (long) (1.5f * Prefs->GetQD3DAtomQuality()),
4251 (long) (Prefs->GetQD3DAtomQuality()));
4252 #endif
4253
4254 glEndList();
4255
4256 gluDeleteQuadric(quad);
4257
4258 // Set the background color to the user's preference.
4259 RGBColor *BackgroundColor = Prefs->GetBackgroundColorLoc();
4260 float red, green, blue;
4261 red = (float) BackgroundColor->red / 65536;
4262 green = (float) BackgroundColor->green / 65536;
4263 blue = (float) BackgroundColor->blue / 65536;
4264 glClearColor(red, green, blue, 1.0f);
4265
4266 // The foreground color is white or black, whichever has better contrast given
4267 // the background color.
4268 long backMagnitude = BackgroundColor->red + BackgroundColor->green + BackgroundColor->blue;
4269 if (backMagnitude > 70000)
4270 fg_color[0] = fg_color[1] = fg_color[2] = 0.0f;
4271 else
4272 fg_color[0] = fg_color[1] = fg_color[2] = 1.0f;
4273
4274 // Setup two lights on either side of the camera.
4275 float fillBrightness = Prefs->GetQD3DFillBrightness();
4276 float PointBrightness = Prefs->GetQD3DPointBrightness();
4277 GLfloat diffuse[4] = {fillBrightness, fillBrightness, fillBrightness, 0.0f};
4278 GLfloat specular[4] = {PointBrightness, PointBrightness, PointBrightness, 0.0f};
4279 GLfloat ambient[4] = {0.2f, 0.2f, 0.2f, 1.0};
4280
4281 // Make sure we're in eye space so that the lights always
4282 // illuminate the visible scene.
4283 glMatrixMode(GL_MODELVIEW);
4284 glLoadIdentity();
4285
4286 glLightfv(GL_LIGHT0, GL_AMBIENT, ambient);
4287 glLightfv(GL_LIGHT0, GL_DIFFUSE, diffuse);
4288 glLightfv(GL_LIGHT0, GL_SPECULAR, specular);
4289 glLightfv(GL_LIGHT0, GL_POSITION, light_pos);
4290 glEnable(GL_LIGHT0);
4291
4292 /* position[0] = -6.0; */
4293 /* ambient[0] = ambient[1] = ambient[2] = 0.0f; */
4294 /* glLightfv(GL_LIGHT1, GL_AMBIENT, ambient); */
4295 /* glLightfv(GL_LIGHT1, GL_DIFFUSE, diffuse); */
4296 /* glLightfv(GL_LIGHT1, GL_SPECULAR, specular); */
4297 /* glLightfv(GL_LIGHT1, GL_POSITION, position); */
4298 /* glEnable(GL_LIGHT1); */
4299
4300 if (GLEW_VERSION_2_0) {
4301 if (depth_fbo) {
4302 glDeleteFramebuffersEXT(1, &(depth_fbo));
4303 depth_fbo = 0;
4304 }
4305
4306 wxString vector_font, bitmap_font;
4307 wxStandardPathsBase & gStdPaths = wxStandardPaths::Get();
4308 #if wxCHECK_VERSION(2, 8, 0)
4309 wxString pathname = gStdPaths.GetResourcesDir();
4310 #else
4311 wxString pathname = gStdPaths.GetDataDir();
4312 #ifdef __WXMAC__
4313 //wxWidgets has a funny idea of where the resources are stored. It locates them as "SharedSupport"
4314 //but xcode is putting them in Resources.
4315 pathname.Remove(pathname.Length() - 13);
4316 pathname += wxT("Resources");
4317 #endif
4318 #endif
4319
4320 std::string vpath, fpath;
4321 if (Prefs->GetShaderMode() == 1) {
4322 vpath = std::string(pathname.ToAscii()) + "/perpixel_dirlight_v.glsl";
4323 fpath = std::string(pathname.ToAscii()) + "/perpixel_dirlight_f.glsl";
4324 shader_program = GetShaderProgramFromFiles(vpath, fpath);
4325 } else if (Prefs->GetShaderMode() == 2) {
4326 vpath = std::string(pathname.ToAscii()) + "/shadows_v.glsl";
4327 fpath = std::string(pathname.ToAscii()) + "/shadows_f.glsl";
4328 shader_program = GetShaderProgramFromFiles(vpath, fpath);
4329
4330 glGenTextures(1, &(depth_tex_id));
4331 glBindTexture(GL_TEXTURE_2D, depth_tex_id);
4332 glTexImage2D(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT16, FBO_SIZE,
4333 FBO_SIZE, 0, GL_DEPTH_COMPONENT, GL_FLOAT, NULL);
4334 glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
4335 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
4336 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
4337 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
4338 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
4339 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
4340 glTexParameteri(GL_TEXTURE_2D, GL_DEPTH_TEXTURE_MODE, GL_INTENSITY);
4341 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_COMPARE_R_TO_TEXTURE);
4342 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_FUNC, GL_LEQUAL);
4343
4344 glGenFramebuffersEXT(1, &depth_fbo);
4345 glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, depth_fbo);
4346 glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT,
4347 GL_DEPTH_ATTACHMENT_EXT, GL_TEXTURE_2D,
4348 depth_tex_id, 0);
4349 glDrawBuffer(GL_NONE);
4350 glReadBuffer(GL_NONE);
4351
4352 GLenum status = glCheckFramebufferStatusEXT(GL_FRAMEBUFFER_EXT);
4353 if (status != GL_FRAMEBUFFER_COMPLETE_EXT) {
4354 std::cout << "Bad status: << " << (int) status << std::endl;
4355 }
4356
4357 glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
4358 }
4359
4360 if (shader_program) {
4361 glUseProgram(shader_program);
4362
4363 CPoint3D lightvec(light_pos[0], light_pos[1], light_pos[2]);
4364 Normalize3D(&lightvec);
4365 GLint light_pos_loc = glGetUniformLocation(shader_program, "light_dir");
4366 if (light_pos_loc >= 0) {
4367 glUniform3f(light_pos_loc, lightvec.x, lightvec.y, lightvec.z);
4368 } else {
4369 std::cerr << "Can't set shader uniforms." << std::endl;
4370 }
4371
4372 GLint uniform_loc = glGetUniformLocation(shader_program, "depth_map");
4373 if (uniform_loc >= 0) {
4374 glUniform1i(uniform_loc, 0);
4375 }
4376
4377 glUseProgram(0);
4378 }
4379 }
4380
4381 }
4382
4383