atomsGL.cpp 20.5 KB
Newer Older
1
2
3
4
5
#include <math.h>

#include "eprintf.h"
#include "TessShaders.h"
#include "UnitCellShaders.h"
6
#include "markerShaders.h"
7
8
9
10
11
#include "atomsGL.h"
#include "atoms.hpp"
#include "ConfigFile.h"
#include "CompileGLShader.h"
#include "polyhedron.h"
12
#include "Grid.h"
13

14
15
16
17
18
19
20
21
int getAtomTimesteps() 
{
	if (fixedAtoms)
		return 1;
	else
		return TIMESTEPS;
}

22
23
24
25
GLenum atomTexture(GLuint t)
{
	GLenum e;
	//rgh: scale atoms here
26
27
28
29
30
31
32
33
	//in google cardboard, this is called again if the program is running, so leave original or atoms get progresivelly smaller!
	float *a=new float[118*4];
	for (int i = 0; i < 118; i++) {
		a[i*4+0]=atomColours[i][0];
		a[i*4+1]=atomColours[i][1];
		a[i*4+2]=atomColours[i][2];
		a[i*4+3]=atomColours[i][3] * atomScaling;
	}
34
35
36
37
38
	glBindTexture(GL_TEXTURE_2D, t); //atom texture
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
39
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, 118, 1, 0, GL_RGBA, GL_FLOAT, a);
40
41
42
43
44

	glBindTexture( GL_TEXTURE_2D, 0 );
	if ((e = glGetError()) != GL_NO_ERROR) {
		eprintf( "opengl error %d, atomTexture\n", e);
	}
45
	delete [] a;
46
47
48
49
50
51
	return e;
}

//WARNING: This should be called after SetupAtoms
//This means that numAtoms now has the cummulative distribution!
//This should be called after the atom texture is prepared, and therefore has the atomscaling pre-multiplied
52
GLenum SetupAtomsNoTess (GLuint **AtomVAO /*[4]*/, GLuint **AtomVertBuffer/*[3]*/, GLuint **AtomIndexBuffer/*[2]*/)
53
54
	//atoms, cloned atoms
	//rgh: FIXME: add AtomVAO[2] for atom trajectories
55
56
57
58
59
{
	//eprintf ("SetupAtomsNoTess 1");
if (!numAtoms)
		return 0;

60
61
62
63
64
if (!solid) {
	eprintf ("SetupAtomsNoTess, error: no solid defined");
	return 0;
}

65
66
67
68
69
70
//eprintf ("SetupAtomsNoTess 2");
	//for now, render an icosahedron
	//http://prideout.net/blog/?p=48 //public domain code
	//xyz nxnynz u=atom type ; 7 floats
	int e;

71
	int totalatoms=numAtoms[getAtomTimesteps() -1];
72
73
	
//eprintf ("SetupAtomsNoTess 2");
74
	*AtomVAO = new GLuint[4]; //atoms, cloned atoms, unused (bonds use Tess atom positions), trajectories
75
	*AtomIndexBuffer= new GLuint[3];//atoms, cloned atoms, bonds
76
	*AtomVertBuffer = new GLuint[3];//atoms, cloned atoms, trajectories
77

78
	glGenVertexArrays(4, *AtomVAO);
79
	glGenBuffers(2, *AtomIndexBuffer);
80
	glGenBuffers(3, *AtomVertBuffer);
81
82
83
84
85
86
87
88
89
//eprintf ("SetupAtomsNoTess 3");
	glBindVertexArray((*AtomVAO)[0]);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, (*AtomIndexBuffer)[0]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[0]);
//eprintf ("SetupAtomsNoTess 4");
	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glEnableVertexAttribArray(2);
	glDisableVertexAttribArray(3);
90
91
	//eprintf ("SetupAtomsNoTess 5, totalatoms=%d, nVerts=%d", totalatoms, solid->nVerts);
	float *tmp = new float[solid->nVerts * 7 * totalatoms];
92
93
	//eprintf ("SetupAtomsNoTess 6");
#ifdef INDICESGL32		
94
	int *tmpi = new int[solid->nFaces*3 * totalatoms];
95
96
97
	//eprintf ("SetupAtomsNoTess 7");
	int *currenti=tmpi;
#else
98
	unsigned short *tmpi = new unsigned short[solid->nFaces*3 * totalatoms];
99
100
101
102
103
104
	//eprintf ("SetupAtomsNoTess 7B");
	unsigned short *currenti=tmpi;
#endif

	float *current=tmp;
	//eprintf ("Before For 1");
105
	for (int p=0;p<getAtomTimesteps() ;p++) {
106
107
		for (int a = 0; a < numAtoms[p]-(p==0?0:numAtoms[p-1]); a++) {
			const int atomNumber = static_cast<int>(atoms[p][4 * a + 3]);
108
			const float radius = atomRadius(atomNumber)*atomScaling;
109
			for (int i = 0; i < solid->nVerts; i++) { //verts
110
				for (int k = 0; k < 3; k++) {
111
					*current++ = solid->Verts[3 * i + k]* radius +atoms[p][4 * a + k]; //pos
112
113
				}
				for (int k = 0; k < 3; k++) {
114
					*current++ = solid->Verts[3 * i + k]; //normal
115
				}
116
				*current++ = static_cast<float>(atomNumber);
117
			} //i
118
119
			for (int i = 0; i < solid->nFaces * 3; i++)
				*currenti++ = solid->Faces[i] + (a+(p==0?0:numAtoms[p-1]))*solid->nVerts;
120
121
		} //a
	} //p
122
		glBufferData(GL_ARRAY_BUFFER, sizeof(float) *totalatoms* 7 * solid->nVerts, tmp,
123
124
125
126
127
128
129
130
131
132
			GL_STATIC_DRAW);
		if ((e = glGetError()) != GL_NO_ERROR)
			eprintf("opengl error %d, glBufferData, l %d\n", e, __LINE__);

		glBufferData(GL_ELEMENT_ARRAY_BUFFER, 
#ifdef INDICESGL32		
	sizeof(int)
#else
	sizeof(unsigned int)
#endif
133
		* totalatoms * 3 * solid->nFaces, tmpi, GL_STATIC_DRAW);
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150


	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 7*sizeof(float), (const void *)0);
	glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(3*sizeof(float)));
	glVertexAttribPointer(2, 1, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(6 * sizeof(float)));

	if (glGetError() != GL_NO_ERROR)
		eprintf("opengl error attrib pointer 0\n");

	//glBindVertexArray(0);
	//glDisableVertexAttribArray(0);
	delete[] tmp;
	delete[] tmpi;
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf("opengl error %d, end of SetupAtoms, l %d\n", e, __LINE__);

	//FIXME TODO: cloned atoms
151
	tmp = new float[solid->nVerts * 7 * numClonedAtoms];
152
153
154
	current=tmp;
	//eprintf ("SetupAtomsNoTess 6");
#ifdef INDICESGL32		
155
	tmpi = new int[solid->nFaces*3 * numClonedAtoms];
156
157
158
	//eprintf ("SetupAtomsNoTess 7");
	currenti=tmpi;
#else
159
	tmpi = new unsigned short[solid->nFaces*3 * numClonedAtoms];
160
161
162
163
164
165
166
	//eprintf ("SetupAtomsNoTess 7B");
	currenti=tmpi;
#endif
	//eprintf ("Before For 2");

	for (int a = 0; a < numClonedAtoms; a++) {
		const int atomNumber = static_cast<int>(clonedAtoms[0][4 * a + 3]);
167
		const float radius = atomRadius(atomNumber)*atomScaling;
168
		for (int i = 0; i < solid->nVerts; i++) { //verts
169
				for (int k = 0; k < 3; k++) {
170
					*current++ = solid->Verts[3 * i + k]* radius +clonedAtoms[0][4 * a + k]; //pos
171
172
				}
				for (int k = 0; k < 3; k++) {
173
					*current++ = solid->Verts[3 * i + k]; //normal
174
				}
175
				*current++ =  static_cast<float>(atomNumber);
176
		} //i
177
178
		for (int i = 0; i < solid->nFaces * 3; i++)
			*currenti++ = solid->Faces[i] + a*solid->nVerts;
179
180
181
182
183
184
185
186
187
188
189
190
191
192
	} //a
	
	//eprintf ("After For 2");


	glBindVertexArray((*AtomVAO)[1]);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, (*AtomIndexBuffer)[1]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[1]);

	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glEnableVertexAttribArray(2);
	glDisableVertexAttribArray(3);

193
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) *numClonedAtoms* 7 * solid->nVerts, tmp,
194
195
196
197
198
199
200
201
202
203
			GL_STATIC_DRAW);
		if ((e = glGetError()) != GL_NO_ERROR)
			eprintf("opengl error %d, glBufferData, l %d\n", e, __LINE__);
	//eprintf ("After bufferdata, array buffer");
		glBufferData(GL_ELEMENT_ARRAY_BUFFER, 
#ifdef INDICESGL32		
	sizeof(int)
#else
	sizeof(unsigned int)
#endif
204
		* numClonedAtoms * 3 * solid->nFaces, tmpi, GL_STATIC_DRAW);
205
206
207
208
209
210
211
212
213
214
	//eprintf ("After bufferdata, element array buffer");
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 7*sizeof(float), (const void *)0);
	glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(3*sizeof(float)));
	glVertexAttribPointer(2, 1, GL_FLOAT, GL_FALSE, 7 * sizeof(float), (const void *)(6 * sizeof(float)));

	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf("opengl error %d, glVertexAttribPointer, l %d\n", e, __LINE__);

	delete[] tmp;
	delete[] tmpi;
215

216
	glBindVertexArray(0);
217
218
219
220
	return e;
} //SetupAtomsNoTess


221
GLenum SetupAtoms(GLuint **AtomVAO /*[4]*/, GLuint **AtomVertBuffer /*[3]*/, GLuint *BondIndices)
222
223
{
	if (!numAtoms)
224
		return glGetError();
225
226
227
228
229
230
231
	//rgh FIXME: put this all in the same vao
	
	//http://prideout.net/blog/?p=48 //public domain code
	//xyz u=atom type ; 4 floats
	int e;

	int totalatoms=0;
232
	for (int i=0;i<getAtomTimesteps() ;i++) {
233
234
235
236
		totalatoms += numAtoms[i];
	}
	eprintf("SetupAtoms: totalatoms=%d", totalatoms);

237
238
	*AtomVAO = new GLuint[4]; //atoms, cloned atoms, bonds, trajectories
	*AtomVertBuffer = new GLuint[3]; //atoms, cloned atoms, trajectories
239

240
241
	glGenVertexArrays(4, *AtomVAO);
	glGenBuffers(3, *AtomVertBuffer);
242
	glGenBuffers(1, BondIndices);
243
244
245
246
247
248
249
250

	glBindVertexArray((*AtomVAO)[0]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[0]);

	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glDisableVertexAttribArray(2);
	glDisableVertexAttribArray(3);
251
252
253
254

	e=glGetError();
	if (e!=GL_NO_ERROR)
		eprintf ("gl error %d, %s %d", e, __FILE__, __LINE__);
255
256
	float *tmp = new float[4 * totalatoms];
	float *current=tmp;
257
258
	
	const int atomlimit=30;
259
	const float bondscaling=0.7f;
260

261
262
	numBonds=new int[getAtomTimesteps() ];
	for (int p=0;p<getAtomTimesteps() ;p++) {
263

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
			for (int a = 0; a < numAtoms[p]; a++) {
				for (int k = 0; k < 4; k++) {
					*current++ = atoms[p][4 * a + k];
				}
			} //a

		if (numAtoms[0]<atomlimit) {
		//eprintf ("searching bonds basic");
		//bonds FIXME quadractic complexity	
				for (int a1=0; a1 < numAtoms[p]; a1++) {
					for (int a2=a1+1; a2 < numAtoms[p]; a2++) {
						float d=0, r;
						for (int k=0;k<3;k++) {
							float dif=atoms[p][4 * a1 + k]-atoms[p][4 * a2 + k];
							d+=dif*dif;
						}
						r=atomRadius(static_cast<int>(atoms[p][4 * a1 + 3]))+
							atomRadius(static_cast<int>(atoms[p][4 * a2 + 3]));
282
						if (d*bondscaling<r*r) {// bond
283
284
285
286
287
288
289
290
291
292
293
							bonds.push_back(a1+(p==0?0:numAtoms[p-1]));
							bonds.push_back(a2+(p==0?0:numAtoms[p-1]));
						}
					}
				}
		} else { //more than 30 atoms, try grid optimization
		//eprintf ("searching bonds grid");

			float m[3];
			float M[3];
			for (int k=0; k<3;k++) {
294
				m[k]=M[k]=atoms[p][k];
295
			}
296
297
			for (int a = 1; a < numAtoms[p]; a++) {
				for (int k=0; k<3;k++) {
298
299
300
301
					if (m[k]>atoms[p][4*a+k])
						m[k]=atoms[p][4*a+k];
					if (M[k]<atoms[p][4*a+k])
						M[k]=atoms[p][4*a+k];
302
303
				}
			}
304
			grid g(m, M, pow(numAtoms[p], 1.0/3), bondscaling);
305
			for (int a = 1; a < numAtoms[p]; a++) 
306
				g.add(atoms[p]+4*a);
307
			for (int a = 0; a < numAtoms[p]; a++) {
308
				std::vector<float*> found=g.find(atoms[p]+4*a);
309
310
311
312
				for (int b=0;b<found.size();b++) {
					//if (found[b] < tmp+4*a) // already got this bound
					//	continue;
					bonds.push_back(a+(p==0?0:numAtoms[p-1]));
313
					bonds.push_back(((found[b]-atoms[p])/4)+(p==0?0:numAtoms[p-1]));
314
315
316
317
				}
			}
		}
		numBonds[p]=bonds.size();
318
319
320
		if (p!=0)
			numAtoms[p]+=numAtoms[p-1];
	} //p
321

322
323
324
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(0));
	glVertexAttribPointer(1, 1, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(3 * sizeof(float)));
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * totalatoms * 4 , tmp,
325
326
327
		GL_STATIC_DRAW);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glBufferData, l %d\n", e, __LINE__);
328

329
	glBindVertexArray(0);
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347

	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, end of SetupAtoms, l %d\n", e, __LINE__);

	if (showTrajectories) {
			//fill the restart buffer
		//use abc for measuring
		float max=0;
		if (has_abc) {
			for (int i=0;i<3;i++)
				for (int j=0;j<3;j++)
				max+=abc[i][j];
			max /=9*2;
		}

		for (unsigned int t=0;t<atomtrajectories.size();t++) {
			atomtrajectoryrestarts.push_back(std::vector<int>());
			atomtrajectoryrestarts[t].push_back(0);
348
			for (int p=1;p<getAtomTimesteps() ;p++) {
349
				int a=atomtrajectories[t];
350
351
352
353
354
				if (has_abc)
					if (fabs(atoms[p][a*4+0]-atoms[p-1][a*4+0])+
						fabs(atoms[p][a*4+1]-atoms[p-1][a*4+1])+
						fabs(atoms[p][a*4+2]-atoms[p-1][a*4+2])>max)
							atomtrajectoryrestarts[t].push_back(p);
355
			}
356
			atomtrajectoryrestarts[t].push_back(getAtomTimesteps() );
357
		}
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
	//need to setup a specific buffer because of GL_MAX_VERTEX_ATTRIB_STRIDE
	//only need xyz, not atom size
	//rgh FIXME: If we use index buffer instead, GPU storage is 1/3 of this
		float *traj = new float[atomtrajectories.size()*TIMESTEPS*3];
		for (unsigned int t = 0; t < atomtrajectories.size(); t++) {
			for (int i=0;i<TIMESTEPS;i++)
				for (int j = 0; j < 3; j++) {
					traj[t*TIMESTEPS * 3 + i * 3 + j] = tmp[i*numAtoms[0]*4+
																+atomtrajectories[t]*4
																+j];
				}
		}
		glBindVertexArray((*AtomVAO)[3]);
		glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[2]);
		glBufferData(GL_ARRAY_BUFFER, sizeof(float) *atomtrajectories.size()*TIMESTEPS * 3, traj,
			GL_STATIC_DRAW);
		glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(float), (const void *)(0));
		glEnableVertexAttribArray(0);
		e = glGetError();
		if ((e = glGetError()) != GL_NO_ERROR)
			eprintf("opengl error %d, creating atom trajectories, l %d\n", e, __LINE__);

		delete[] traj;
381
382
	}
	delete[] tmp;
383
384
385
386
387
388
389
390
391
392
393
394
	//bonds
	glBindVertexArray((*AtomVAO)[2]);
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[0]);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, *BondIndices);
	glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(int)*bonds.size(), bonds.data(), GL_STATIC_DRAW);
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(0));
	glEnableVertexAttribArray(0);
	glBindVertexArray(0);

	e=glGetError();
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, creating chemical bonds, l %d\n", e, __LINE__);
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414

	//now clones
	if (basisvectorreps ||!clonedAtoms) //do not replicate
		return e;


	glBindVertexArray((*AtomVAO)[1]); //rgh FIXME, only works for TIMESTEPS=1
	glBindBuffer(GL_ARRAY_BUFFER, (*AtomVertBuffer)[1]);
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * clonedAtoms[0].size(), clonedAtoms[0].data(),
			GL_STATIC_DRAW);
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(0));
	glVertexAttribPointer(1, 1, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (const void *)(3 * sizeof(float)));
	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, end of Setup cloned Atoms, l %d\n", e, __LINE__);

	//rgh: we will need these again if we don't have tesselation
	//delete[] clonedAtoms;
	//clonedAtoms=0;
415
	glBindVertexArray(0);
416
417
418
	return e;
}

419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
GLenum SetupMarker(GLuint *MarkerVAO, GLuint *MarkerVertBuffer)
{
	if (!markers)
		return glGetError();
	GLenum e;
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, begin of SetupMarker\n", e, __LINE__);

	glGenVertexArrays(1, MarkerVAO);
	glGenBuffers(1, MarkerVertBuffer);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glGenBuffers, l %d\n", e, __LINE__);

	glBindVertexArray(*MarkerVAO);
	glBindBuffer(GL_ARRAY_BUFFER, *MarkerVertBuffer);

	glEnableVertexAttribArray(0);
	glEnableVertexAttribArray(1);
	glDisableVertexAttribArray(2);
	glDisableVertexAttribArray(3);

	const float size=atomRadius(0)*atomScaling*markerscaling;
	float *tmp = new float [8*TIMESTEPS];
	for (int i=0;i<TIMESTEPS;i++) {
		for (int j=0;j<3;j++) { //center [3]
			tmp[i*8+j]=markers[i][j];
		}
		tmp[i*8+3]=0.8*size; //size [1]
		for (int j=0;j<4;j++) {//colour[4]
			tmp[i*8+4+j]=markercolours[i][j];
		}
	}
	
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * TIMESTEPS*8 , tmp,
			GL_STATIC_DRAW);
	glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 8 * sizeof(float), (const void *)(0));
	glVertexAttribPointer(1, 4, GL_FLOAT, GL_FALSE, 8 * sizeof(float), (const void *)(4*sizeof(float)));
	glBindVertexArray(0);
	return glGetError();
}

460
461
GLenum SetupUnitCell(GLuint *UnitCellVAO, GLuint *UnitCellVertBuffer, GLuint *UnitCellIndexBuffer)
{
462
	//add here both unit cell and supercell
463
	GLenum e;
464
465
466
467
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, begin of SetupUnitCell\n", e, __LINE__);
	if (!has_abc)
		return e;
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
	glGenVertexArrays(1, UnitCellVAO);
	glGenBuffers(1, UnitCellVertBuffer);
	glGenBuffers(1, UnitCellIndexBuffer);
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glGenBuffers, l %d\n", e, __LINE__);

	glBindVertexArray(*UnitCellVAO);
	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, *UnitCellIndexBuffer);
	glBindBuffer(GL_ARRAY_BUFFER, *UnitCellVertBuffer);

	glEnableVertexAttribArray(0);
	glDisableVertexAttribArray(1);
	glDisableVertexAttribArray(2);
	glDisableVertexAttribArray(3);

483
	float *tmp = new float[3*8*2];
484
	//0, a, b, c, a+b+c, b+c, a+c, a+b
485
	for (int i=0;i<3;i++) { //unit cell
486
487
488
489
490
491
492
493
		tmp[0+i]=0;
		for (int j=0;j<3;j++)
			tmp[3*(j+1)+i]=abc[j][i];
		tmp[3*4+i]=abc[0][i]+abc[1][i]+abc[2][i];
		tmp[3*5+i]=			abc[1][i]+abc[2][i];
		tmp[3*6+i]=abc[0][i]+		abc[2][i];
		tmp[3*7+i]=abc[0][i]+abc[1][i];
	}
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
	float displ[3]={0,0,0};
	if (translations && ISOS) 
		for (int j=0;j<3;j++)
			for (int i=0;i<3;i++)
				displ[i]+=-translations[0][i]*abc[j][i];
	for (int i=0;i<3;i++) { //rgh fixme, add displacement here as well
		tmp[3*8+i]=displ[i];
		for (int j=0;j<3;j++)
			tmp[3*(j+8+1)+i]=abc[j][i]*supercell[j]+displ[i];
		tmp[3*12+i]=abc[0][i]*supercell[0]+abc[1][i]*supercell[1]+abc[2][i]*supercell[2]+displ[i];
		tmp[3*13+i]=			abc[1][i]*supercell[1]+abc[2][i]*supercell[2]+displ[i];
		tmp[3*14+i]=abc[0][i]*supercell[0]+		abc[2][i]*supercell[2]+displ[i];
		tmp[3*15+i]=abc[0][i]*supercell[0]+abc[1][i]*supercell[1]+displ[i];
	}
	int tmpi[12*2*2]={ //lines, unit cell, 
509
510
511
512
513
514
515
516
517
518
519
		0,1, 
		1,6,
		6,3,
		3,0,
		2,7,
		7,4,
		4,5,
		5,2,
		0,2,
		1,7,
		6,4,
520
521
522
523
524
525
526
527
528
529
530
531
532
		3,5, // supercell
		0+8,1+8, 
		1+8,6+8,
		6+8,3+8,
		3+8,0+8,
		2+8,7+8,
		7+8,4+8,
		4+8,5+8,
		5+8,2+8,
		0+8,2+8,
		1+8,7+8,
		6+8,4+8,
		3+8,5+8
533
	};
534
	glBufferData(GL_ARRAY_BUFFER, sizeof(float) * 3*8*2 , tmp,
535
536
			GL_STATIC_DRAW);
	if ((e = glGetError()) != GL_NO_ERROR)
537
		eprintf( "opengl error %d, glBufferData vertex, l %d\n", e, __LINE__);
538
539
	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(float), (const void *)(0));
	glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(tmpi), tmpi, GL_STATIC_DRAW);
540
541
	if ((e = glGetError()) != GL_NO_ERROR)
		eprintf( "opengl error %d, glBufferData index, l %d\n", e, __LINE__);
542
	glBindVertexArray(0);
543
544
545
546
	return e;
}


547
548
549
bool PrepareUnitCellAtomShader (GLuint *AtomP, GLuint *cellP, GLuint *MarkerP, 
								GLint *AtomMatrixLocation, GLint *UnitCellMatrixLocation,  GLint *UnitCellColourLocation,
								GLint *MarkerMatrixLocation){
550
551
552
553
554
555
	if (!PrepareAtomShader(AtomP, AtomMatrixLocation))
		return false;

	if (!PrepareUnitCellShader(cellP, UnitCellMatrixLocation, UnitCellColourLocation))
		return false;

556
557
558
	if (!PrepareMarkerShader(MarkerP, MarkerMatrixLocation))
		return false;

559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
	return true;
}

bool PrepareAtomShader (GLuint *AtomP, GLint *AtomMatrixLocation){
		//https://www.gamedev.net/topic/591110-geometry-shader-point-sprites-to-spheres/
	//no rotation, only translations means we can do directional lighting in the shader.
	//FIXME
	//http://stackoverflow.com/questions/40101023/flat-shading-in-webgl
	*AtomP = CompileGLShader(
		AtomShaders[SHADERNAME],
		AtomShaders[SHADERVERTEX],
		AtomShaders[SHADERFRAGMENT],
		AtomShaders[SHADERTESSEVAL]
		);
	*AtomMatrixLocation=glGetUniformLocation(*AtomP, "matrix");
	if( *AtomMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in atom shader\n" );
		return false;
	}
	return true;
}

582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
bool PrepareMarkerShader (GLuint *MP, GLint *MMatrixLocation){
		//https://www.gamedev.net/topic/591110-geometry-shader-point-sprites-to-spheres/
	//no rotation, only translations means we can do directional lighting in the shader.
	//FIXME
	//http://stackoverflow.com/questions/40101023/flat-shading-in-webgl
	*MP = CompileGLShader(
		MarkerShaders[SHADERNAME],
		MarkerShaders[SHADERVERTEX],
		MarkerShaders[SHADERFRAGMENT],
		MarkerShaders[SHADERTESSEVAL]
		);
	*MMatrixLocation=glGetUniformLocation(*MP, "matrix");
	if( *MMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in atom shader\n" );
		return false;
	}
	return true;
}

602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
bool PrepareAtomShaderNoTess (GLuint *AtomP, GLint *AtomMatrixLocation){
		//https://www.gamedev.net/topic/591110-geometry-shader-point-sprites-to-spheres/
	//no rotation, only translations means we can do directional lighting in the shader.
	//FIXME
	//http://stackoverflow.com/questions/40101023/flat-shading-in-webgl
	*AtomP = CompileGLShader(
		AtomShadersNoTess[SHADERNAME],
		AtomShadersNoTess[SHADERVERTEX],
		AtomShadersNoTess[SHADERFRAGMENT],
		AtomShadersNoTess[SHADERTESSEVAL]
		);
	*AtomMatrixLocation=glGetUniformLocation(*AtomP, "matrix");
	if( *AtomMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in atom shader no tess\n" );
		return false;
	}
	return true;
}


bool PrepareUnitCellShader (GLuint *cellP, GLint *UnitCellMatrixLocation,  GLint *UnitCellColourLocation){
	*cellP= CompileGLShader(
		UnitCellShaders[SHADERNAME],
		UnitCellShaders[SHADERVERTEX],
		UnitCellShaders[SHADERFRAGMENT],
		UnitCellShaders[SHADERTESSEVAL]
		);
	*UnitCellMatrixLocation=glGetUniformLocation(*cellP, "matrix");
	if( *UnitCellMatrixLocation == -1 )
	{
		eprintf( "Unable to find matrix uniform in UnitCell shader\n" );
		return false;
	}
	*UnitCellColourLocation=glGetUniformLocation(*cellP, "color");
	if( *UnitCellColourLocation == -1 )
	{
		eprintf( "Unable to find color uniform in UnitCell shader\n" );
		return false;
	}
	return true;
}


/**p: input, f: output*/
void GetDisplacement(int p[3], float f[3])
{
float delta[3][3];
for (int ss=0;ss<3;ss++)
	for (int i=0;i<3;i++)
		delta[ss][i]=static_cast<float>(p[ss])*abc[ss][i];

for (int i=0;i<3;i++)
	f[i]=0;

for (int ss=0;ss<3;ss++)
	for (int i=0;i<3;i++)
		f[i]+=delta[ss][i];
660
}