1 /*
2 * nsvc.cpp
3 *
4 * Copyright 2012 Martin Robinson
5 *
6 * This file is part of Smoldyn.
7 *
8 * Smoldyn is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU Lesser General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * Smoldyn is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public License
19 * along with Smoldyn. If not, see <http://www.gnu.org/licenses/>.
20 *
21 * Created on: Apr 16, 2013
22 * Author: mrobins
23 */
24
25 #define NSVC_CPP
26 #include "nsvc.h"
27 #include "smoldynfuncs.h"
28 #include <sstream>
29 #include <set>
30 #include <numeric>
31
nsv_init()32 void nsv_init() {
33 Kairos::init(0,NULL);
34 }
35
nsv_new(double * min,double * max,double * dx,int n)36 NextSubvolumeMethod* nsv_new(double* min, double* max, double* dx, int n) {
37 using namespace Kairos;
38 Vect3d min_vect(0.0,0.0,0.0);
39 Vect3d max_vect(1.0,1.0,1.0);
40 Vect3d dx_vect(1.0,1.0,1.0);
41 for (int i = 0; i < n; ++i) {
42 min_vect[i] = min[i];
43 max_vect[i] = max[i];
44 dx_vect[i] = dx[i];
45 }
46 StructuredGrid *grid = new StructuredGrid(min_vect,max_vect,dx_vect);
47 NextSubvolumeMethod* nsv = new NextSubvolumeMethod(*grid);
48
49 return nsv;
50 }
51
nsv_delete(NextSubvolumeMethod * nsv)52 void nsv_delete(NextSubvolumeMethod* nsv) {
53 if (!nsv) return;
54 delete &nsv->get_grid();
55 delete nsv;
56 }
57
nsv_print(NextSubvolumeMethod * nsv,char ** bufferptr)58 void nsv_print(NextSubvolumeMethod* nsv, char** bufferptr) {
59 int length;
60 std::ostringstream ss;
61 char *buffer;
62
63 ss << std::endl << *nsv << std::endl;
64 length=ss.str().length();
65 buffer=(char*) calloc(length,sizeof(char));
66 if(!buffer) return;
67 ss.str().copy(buffer,length,0);
68 buffer[length-1] = '\0';
69 *bufferptr=buffer;
70 return;
71 }
72
73
nsv_add_interface(NextSubvolumeMethod * nsv,int id,double dt,double * start,double * end,double * norm,int dim)74 void nsv_add_interface(NextSubvolumeMethod* nsv,int id,double dt, double *start,double *end,double *norm,int dim) {
75 using namespace Kairos;
76 Vect3d min(0.0,0.0,0.0);
77 Vect3d max(1.0,1.0,1.0);
78 for (int i = 0; i < dim; ++i) {
79 min[i] = start[i];
80 max[i] = end[i];
81 }
82 if ((norm[0] == 1) || (norm[0] == -1)) {
83 xrect interface(min,max,norm[0]);
84 nsv->set_interface(interface,id,dt,false);
85
86 } else if ((norm[1] == 1) || (norm[1] == -1)) {
87 yrect interface(min,max,norm[1]);
88 nsv->set_interface(interface,id,dt,false);
89
90 } else if ((norm[2] == 1) || (norm[2] == -1)) {
91 zrect interface(min,max,norm[2]);
92 nsv->set_interface(interface,id,dt,false);
93
94 }
95 }
96
97
nsv_add_species(NextSubvolumeMethod * nsv,int id,double D,char * btype,int dim)98 void nsv_add_species(NextSubvolumeMethod* nsv,int id,double D,char *btype,int dim) {
99 using namespace Kairos;
100
101 //TODO: assumes that species has not been added before
102 Species *ns = new Species(id,D,&nsv->get_grid());
103 nsv->add_diffusion(*ns);
104
105 if (dim > 0) {
106 xplane xlow(nsv->get_grid().get_low()[0],1);
107 xplane xhigh(nsv->get_grid().get_high()[0],-1);
108 if (btype[0]=='p') {
109
110 const double hx = nsv->get_grid().get_cell_size()[0];
111 nsv->add_diffusion_between(*ns,ns->D/pow(hx,2),xlow,xhigh);
112 nsv->add_diffusion_between(*ns,ns->D/pow(hx,2),xhigh,xlow);
113 }
114 if (dim > 1) {
115 yplane ylow(nsv->get_grid().get_low()[1],1);
116 yplane yhigh(nsv->get_grid().get_high()[1],-1);
117 if (btype[1]=='p') {
118
119 const double hy = nsv->get_grid().get_cell_size()[1];
120 nsv->add_diffusion_between(*ns,ns->D/pow(hy,2),ylow,yhigh);
121 nsv->add_diffusion_between(*ns,ns->D/pow(hy,2),yhigh,ylow);
122 }
123 if (dim > 2) {
124 zplane zlow(nsv->get_grid().get_low()[2],1);
125 zplane zhigh(nsv->get_grid().get_high()[2],-1);
126 if (btype[2]=='p') {
127
128 const double hz = nsv->get_grid().get_cell_size()[2];
129 nsv->add_diffusion_between(*ns,ns->D/pow(hz,2),zlow,zhigh);
130 nsv->add_diffusion_between(*ns,ns->D/pow(hz,2),zhigh,zlow);
131 }
132 }
133 }
134 }
135 }
136
137 class SmoldynSurface {
138 public:
SmoldynSurface(surfaceptr s,enum PanelFace face)139 SmoldynSurface(surfaceptr s, enum PanelFace face):s(s),face(face) {}
lineXsurface(const Kairos::Vect3d & p1,const Kairos::Vect3d & p2) const140 bool lineXsurface(const Kairos::Vect3d& p1, const Kairos::Vect3d& p2) const {
141 if (face==PFnone) return false;
142 double crsspt[3];
143 double p1t[3],p2t[3];
144 enum PanelFace face1,face2;
145 p1t[0] = p1[0];
146 p1t[1] = p1[1];
147 p1t[2] = p1[2];
148 p2t[0] = p2[0];
149 p2t[1] = p2[1];
150 p2t[2] = p2[2];
151 for (int i = 0; i < PSMAX; ++i) {
152 for (int j = 0; j < s->npanel[i]; ++j) {
153 if (lineXpanel(p1t,p2t,s->panels[i][j],s->srfss->sim->dim,crsspt,&face1,&face2,NULL,NULL,NULL,0)==1) {
154 if ((face1!=face2) && ((face==PFboth)||(face1==face))) return true;
155 }
156 }
157 }
158 return false;
159 }
160
161
162 private:
163 surfaceptr s;
164 enum PanelFace face;
165 };
166
167
168 class SmoldynCompartment {
169 public:
SmoldynCompartment(compartptr c)170 SmoldynCompartment(compartptr c):c(c) {}
is_in(const Kairos::Vect3d & test_point) const171 bool is_in(const Kairos::Vect3d& test_point) const {
172 double pt[3];
173 pt[0] = test_point[0];
174 pt[1] = test_point[1];
175 pt[2] = test_point[2];
176 return posincompart(c->cmptss->sim,pt,c,0)==1;
177 }
178 private:
179 compartptr c;
180 };
181
nsv_add_surface(NextSubvolumeMethod * nsv,surfacestruct * surface)182 extern void nsv_add_surface(NextSubvolumeMethod* nsv,surfacestruct* surface) {
183 using namespace Kairos;
184 const int ns = nsv->get_diffusing_species().size();
185 //std::cout << "adding surface "<<surface->sname<<std::endl;
186 for (int i = 0; i < ns; ++i) {
187 enum PanelFace face;
188 int id = nsv->get_diffusing_species()[i]->id;
189 if ((surface->action[id][0][PFfront]==SAreflect) &&
190 (surface->action[id][0][PFback]==SAreflect)) face = PFboth;
191 else if (surface->action[id][0][PFfront]==SAreflect) face = PFfront;
192 else if (surface->action[id][0][PFback]==SAreflect) face = PFback;
193 else face = PFnone;
194 if (face != PFnone) {
195 //std::cout << "finished reflective surface with face="<<face<<std::endl;
196 SmoldynSurface s(surface,face);
197 nsv->scale_diffusion_across(*(nsv->get_diffusing_species()[i]),s,0);
198 }
199
200 if ((surface->action[id][0][PFfront]==SAabsorb) &&
201 (surface->action[id][0][PFback]==SAabsorb)) face = PFboth;
202 else if (surface->action[id][0][PFfront]==SAabsorb) face = PFfront;
203 else if (surface->action[id][0][PFback]==SAabsorb) face = PFback;
204 else face = PFnone;
205 if (face != PFnone) {
206 //std::cout << "finished absorption surface with face="<<face<<std::endl;
207 SmoldynSurface s(surface,face);
208 nsv->absorption_across(*(nsv->get_diffusing_species()[i]),s,1);
209 }
210 }
211 //std::cout << "finished adding surface "<<surface->sname<<std::endl;
212
213 }
214
215
nsv_add_reaction(NextSubvolumeMethod * nsv,rxnstruct * reaction)216 void nsv_add_reaction(NextSubvolumeMethod* nsv,rxnstruct *reaction) {
217 const double rate = reaction->rate;
218 const int nreactants = reaction->rxnss->order;
219 const int *reactant_ids = reaction->rctident;
220 const int nproducts = reaction->nprod;
221 const int *product_ids = reaction->prdident;
222
223
224 using namespace Kairos;
225
226
227 ReactionSide lhs;
228 for (int i = 0; i < nreactants; ++i) {
229 Species* s = nsv->get_species(reactant_ids[i]);
230 bool found = false;
231 for (std::vector<ReactionComponent>::iterator rc=lhs.begin();rc!=lhs.end();rc++) {
232 if (rc->species==s) {
233 rc->multiplier++;
234 found = true;
235 }
236 }
237 if (!found) {
238 lhs = lhs + *(nsv->get_species(reactant_ids[i]));
239 }
240 }
241 ReactionSide rhs;
242 for (int i = 0; i < nproducts; ++i) {
243 Species* s = nsv->get_species(product_ids[i]);
244
245 bool to_particle = false;
246 if (reaction->prdrep != NULL) to_particle = (reaction->prdrep[i] == SRparticle);
247
248 bool found = false;
249 for (std::vector<ReactionComponent>::iterator rc=rhs.begin();rc!=rhs.end();rc++) {
250 bool rc_to_particle = (rc->compartment_index == -1);
251 if ((rc->species==s) && (to_particle == rc_to_particle)) {
252 rc->multiplier++;
253 found = true;
254 }
255 }
256 if (!found) {
257 rhs = rhs + *(nsv->get_species(product_ids[i]));
258 if (to_particle) {
259 // std::cout << "adding a t0_paricle reaction "<<rhs<<std::endl;
260 (rhs.end()-1)->compartment_index = -1;
261 } else {
262 // std::cout << "adding a normal reaction"<<std::endl;
263 (rhs.end()-1)->compartment_index = 1;
264 }
265 }
266 }
267
268 if (reaction->srf) {
269 enum PanelFace face = PFboth;
270 SmoldynSurface surface(reaction->srf,face);
271 nsv->add_reaction_on(rate, lhs >> rhs, surface);
272 } else if (reaction->cmpt) {
273 SmoldynCompartment compartment(reaction->cmpt);
274 nsv->add_reaction_in(rate, lhs >> rhs, compartment);
275 } else {
276 nsv->add_reaction(rate, lhs >> rhs);
277 }
278 }
279
280
281 //void nsv_add_reaction(NextSubvolumeMethod* nsv,double rate,
282 // int nreactants,int *reactant_ids,
283 // int nproducts,int *product_ids) {
284 //
285 // //TODO: assumes that species has already been added via nsv_add_species
286 // using namespace Kairos;
287 //
288 //
289 // ReactionSide lhs;
290 // for (int i = 0; i < nreactants; ++i) {
291 // lhs = lhs + *(nsv->get_species(reactant_ids[i]));
292 // }
293 // ReactionSide rhs;
294 // for (int i = 0; i < nproducts; ++i) {
295 // rhs = rhs + *(nsv->get_species(product_ids[i]));
296 // }
297 // nsv->add_reaction(rate, lhs >> rhs);
298 //}
299
nsv_integrate(NextSubvolumeMethod * nsv,double dt,portstruct * port,latticestruct * lattice)300 void nsv_integrate(NextSubvolumeMethod* nsv,double dt, portstruct *port, latticestruct *lattice) {
301 using namespace Kairos;
302 //std::cout << "running lattice dt"<<std::endl;
303 nsv->integrate(dt);
304 const int ns = nsv->get_diffusing_species().size();
305 int n = 0;
306 for (int i = 0; i < ns; ++i) {
307 n += nsv->get_diffusing_species()[i]->particles.size();
308 }
309
310 int *species = new int[n];
311 double **positions = new double*[n];
312 double **positionsx = new double*[n];
313
314 //const double lattice_lengthscale = 0.0001*std::min(lattice->dx[0],std::min(lattice->dx[1],lattice->dx[2]));
315 if (!port) return;
316
317 simptr sim = port->portss->sim;
318 n += sim->mols->nl[port->llport];
319 double *crsspt = new double[sim->dim];
320 enum PanelFace face1,face2;
321
322 /*
323 * look in port for new particles
324 */
325 std::set<int> dirty_indicies;
326 int nout = 0;
327 n = sim->mols->nl[port->llport];
328 for (int i = 0; i < n; ++i) {
329 moleculeptr m = sim->mols->live[port->llport][i];
330
331 // //process other surface interactions
332 // int er = checksurfaces1mol(sim,m);
333 // if (er != 0) simLog(NULL,11,"ERROR: failure in nsv_integrate (while processing surface interactions)\n");
334 //
335 // //check if already gone back through port
336 //
337 // bool in = true;
338 // for (int j = 0; j < port->srf->npanel[PSrect]; ++j) {
339 // if (lineXpanel(m->posx,m->pos,port->srf->panels[PSrect][j],sim->dim,crsspt,&face1,&face2,NULL,NULL,NULL,0)) {
340 // if ((face1!=face2) && (face2==PFfront)) {
341 // in = false;
342 // }
343 // }
344 // }
345 //
346 // if (!in) {
347 // //std::cout <<" particle going back through port!!!"<<std::endl;
348 // species[nout] = m->ident;
349 // positions[nout] = m->pos;
350 // positionsx[nout] = m->via;
351 // nout++;
352 // continue;
353 // }
354
355
356 // if outside lattice domain raise error
357 Vect3d newr(0.5,0.5,0.5);
358 for (int d = 0; d < sim->dim; ++d) {
359 newr[d] = m->via[d];
360 double low = nsv->get_grid().get_low()[d];
361 double high = nsv->get_grid().get_high()[d];
362 if (newr[d] < low) {
363 std::cout << "d = "<<d<<" via = "<<m->via[d]<<" pos = "<<m->pos[d]<<" posx = "<<m->posx[d]<<" newr = "<<newr[d]<<std::endl;
364 simLog(NULL,11,"ERROR: particle unexpectedly outside lattice domain\n");
365 } else if (newr[d] > high) {
366 std::cout << "d = "<<d<<" via = "<<m->via[d]<<" pos = "<<m->pos[d]<<" posx = "<<m->posx[d]<<" newr = "<<newr[d]<<std::endl;
367 simLog(NULL,11,"ERROR: particle unexpectedly outside lattice domain\n");
368 }
369 }
370
371 const int ci = nsv->get_grid().get_cell_index(newr);
372 Vect3d cc = nsv->get_grid().get_cell_centre(ci);
373
374 bool throw_back = false;
375
376 for (int i = 0; i < lattice->nsurfaces; ++i) {
377 if ((lattice->surfacelist[i]->action[m->ident][0][PFfront] != SAreflect) && (lattice->surfacelist[i]->action[m->ident][0][PFback] != SAreflect)) {
378 continue;
379 }
380 for (int j = 0; j < PSMAX; ++j) {
381 for (int k = 0; k < lattice->surfacelist[i]->npanel[j]; ++k) {
382 if (lineXpanel(m->posx,cc.data(),lattice->surfacelist[i]->panels[j][k],port->portss->sim->dim,crsspt,&face1,&face2,NULL,NULL,NULL,0)==1) {
383 //if particle crossed surface when adding to lattice, throw back through port
384 if ((face1!=face2) && (lattice->surfacelist[i]->action[m->ident][0][face1]==SAreflect)) throw_back = true;
385
386 }
387 }
388 }
389 }
390
391
392 if (throw_back) {
393 //std::cout << "throwing back particle at posx = ("<<m->posx[0]<<','<<m->posx[1]<<','<<m->posx[2]<<") and pos = ("<<m->pos[0]<<','<<m->pos[1]<<','<<m->pos[2]<<")"<<std::endl;
394 species[nout] = m->ident;
395 positions[nout] = m->posx;
396 positionsx[nout] = m->posx;
397 nout++;
398 } else {
399 nsv->get_species(m->ident)->copy_numbers[ci]++;
400 dirty_indicies.insert(ci);
401 }
402 }
403 delete[] crsspt;
404
405
406 //delete particles in port
407 int er = portgetmols(sim,port,-1,MSall,1);
408 if (er != n) simLog(NULL,11,"ERROR: failure in nsv_integrate (while deleting particles from port)\n");
409
410
411 /*
412 * put new particles in port
413 */
414 for (int i = 0; i < ns; ++i) {
415 Species *s = nsv->get_diffusing_species()[i];
416 const int np = s->particles.size();
417 double *crsspt = new double[sim->dim];
418 enum PanelFace face1,face2;
419 for (int j = 0; j < np; ++j) {
420 //check if already gone back through port
421 bool in = true;
422 if (port->srf) {
423 for (int k = 0; k < port->srf->npanel[PSrect]; ++k) {
424 if (lineXpanel(s->particlesx[j].data(),s->particles[j].data(),port->srf->panels[PSrect][k],sim->dim,crsspt,&face1,&face2,NULL,NULL,NULL,0)) {
425 if ((face1!=face2) && (port->srf->action[s->id][0][face1]==SAport)) {
426 in = false;
427 }
428 }
429 }
430 }
431
432 if (in) {
433 //std::cout <<" particle going back through port!!!"<<std::endl;
434 species[nout] = s->id;
435 positions[nout] = s->particles[j].data();
436 // if (positions[nout][0] <= -10) std::cout << "x <= -10: "<<positions[nout]<<std::endl;
437 // if (positions[nout][0] >= 10) std::cout << "x >= -10: "<<positions[nout]<<std::endl;
438 // if (positions[nout][1] <= -10) std::cout << "x <= -10: "<<positions[nout]<<std::endl;
439 // if (positions[nout][1] >= 10) std::cout << "x >= -10: "<<positions[nout]<<std::endl;
440 positionsx[nout] = s->particlesx[j].data();
441 nout++;
442 } else {
443 //add to lattice
444 //std::cout << "adding particle of species "<<s->id<<" back into lattice at position "<<s->particles[j]<<std::endl;
445 const int ci = nsv->get_grid().get_cell_index(s->particles[j]);
446 //Vect3d cc = nsv->get_grid().get_cell_centre(ci);
447 nsv->get_species(s->id)->copy_numbers[ci]++;
448 dirty_indicies.insert(ci);
449 }
450
451 }
452 s->particles.clear();
453 s->particlesx.clear();
454 }
455
456 for (std::set<int>::iterator i=dirty_indicies.begin();i!=dirty_indicies.end();i++) {
457 nsv->recalc_priority(*i);
458 }
459
460 //put particles back in port if needed
461 if (nout > 0) {
462 er = portputmols(sim,port,nout,species[0],species,positions,positionsx);
463 if (er != 0) simLog(NULL,11,"ERROR: failure in nsv_integrate (while putting particles in the port)\n");
464 }
465
466 delete[] species;
467 delete[] positions;
468 }
469
nsv_get_grid(NextSubvolumeMethod * nsv)470 vtkUnstructuredGrid* nsv_get_grid(NextSubvolumeMethod* nsv) {
471 using namespace Kairos;
472 (void)nsv;
473 #if defined(HAVE_VTK)
474 return get_vtk_grid(nsv->get_grid(),nsv->get_diffusing_species());
475 #endif
476 return NULL;
477 }
478
nsv_molcountspace(NextSubvolumeMethod * nsv,int id,double * low,double * high,int dim,int nbins,int axis,int * ret_array)479 void nsv_molcountspace(NextSubvolumeMethod* nsv,int id, double *low, double *high, int dim, int nbins, int axis, int *ret_array) {
480 using namespace Kairos;
481
482 Vect3d vlow(0,0,0);
483 Vect3d vhigh(1,1,1);
484 Vect3d grid_size(1,1,1);
485
486 for (int i = 0; i < dim; ++i) {
487 vlow[i] = low[i];
488 vhigh[i] = high[i];
489 grid_size[i] = high[i] - low[i];
490 }
491 if (nbins > 1) {
492 grid_size[axis] = (high[axis] - low[axis])/nbins;
493 }
494
495 StructuredGrid grid(vlow,vhigh,grid_size);
496
497 std::vector<double> concentration;
498
499 nsv->get_species(id)->get_concentration(grid,concentration);
500
501 for (int i = 0; i < nbins; ++i) {
502 ret_array[i] = concentration[i]*grid.get_cell_volume(i);
503 }
504 }
505
506 /*
507 * returns the concentration of species id at point
508 */
nsv_concentration_point(NextSubvolumeMethod * nsv,int id,double * point,int dim)509 double nsv_concentration_point(NextSubvolumeMethod* nsv,int id, double *point, int dim) {
510 using namespace Kairos;
511 Vect3d vpoint(0,0,0);
512 for (int i = 0; i < dim; ++i) {
513 vpoint[i] = point[i];
514 }
515 const Species* s = nsv->get_species(id);
516 const int index = nsv->get_grid().get_cell_index(vpoint);
517 return s->copy_numbers[index]/nsv->get_grid().get_cell_volume(index);
518 }
519
nsv_molcount(NextSubvolumeMethod * nsv,int * ret_array)520 void nsv_molcount(NextSubvolumeMethod* nsv, int *ret_array) {
521 using namespace Kairos;
522
523 std::vector<Species*> species = nsv->get_diffusing_species();
524 for (unsigned int i = 0; i < species.size(); ++i) {
525 const int n = std::accumulate(species[i]->copy_numbers.begin(),species[i]->copy_numbers.end(),0) +
526 species[i]->particles.size();
527 ret_array[species[i]->id] = n;
528 }
529 }
530
531 /*
532 * kills one molecule of species id at point
533 */
nsv_kill_molecule(NextSubvolumeMethod * nsv,int id,double * point,int dim)534 void nsv_kill_molecule(NextSubvolumeMethod* nsv,int id, double *point, int dim) {
535 using namespace Kairos;
536 Vect3d vpoint(0,0,0);
537 for (int i = 0; i < dim; ++i) {
538 vpoint[i] = point[i];
539 }
540 Species* s = nsv->get_species(id);
541 const int index = nsv->get_grid().get_cell_index(vpoint);
542 s->copy_numbers[index]--;
543 if (s->copy_numbers[index] < 0) {
544 simLog(NULL,11,"ERROR: lattice species became less than zero (in nsv_kill_molecule)\n");
545 }
546 nsv->recalc_priority(index);
547 }
548
549
nsv_add_mol(NextSubvolumeMethod * nsv,int id,double * pos,int dim)550 void nsv_add_mol(NextSubvolumeMethod* nsv,int id, double* pos, int dim) {
551 using namespace Kairos;
552 Vect3d newr(0.5,0.5,0.5);
553 // if outside lattice domain raise error
554 for (int d = 0; d < dim; ++d) {
555 double low = nsv->get_grid().get_low()[d];
556 double high = nsv->get_grid().get_high()[d];
557 if (pos[d] < low) {
558 simLog(NULL,11,"ERROR: particle unexpectedly outside lattice domain\n");
559 } else if (pos[d] > high) {
560 simLog(NULL,11,"ERROR: particle unexpectedly outside lattice domain\n");
561 } else {
562 newr[d] = pos[d];
563 }
564 }
565 //nsv->get_species(id)->particles.push_back(newr);
566
567 const int ci = nsv->get_grid().get_cell_index(newr);
568 nsv->get_species(id)->copy_numbers[ci]++;
569 nsv->recalc_priority(ci);
570 }
571
572 /*
573 * returns vector of copy numbers for each lattice site and positions of lattice sites
574 *
575 * return value = number of lattice sites returned for species "id"
576 * copy_numbers = array of ints representing copy number of species "id" at each lattice site
577 * equals NULL if species does not exist
578 * positions = array of doubles representing positions of each site. In lattice-site major order.
579 * eg. x-coord of lattice site n is positions[3*n+0]
580 */
nsv_get_species_copy_numbers(NextSubvolumeMethod * nsv,int id,const int ** copy_numbers,const double ** positions)581 int nsv_get_species_copy_numbers(NextSubvolumeMethod* nsv, int id,const int **copy_numbers, const double** positions) {
582 using namespace Kairos;
583 int num_lattice_sites;
584
585 Species* s = nsv->get_species(id);
586 if (s==NULL) {
587 *copy_numbers = NULL;
588 num_lattice_sites = 0;
589 } else {
590 *copy_numbers = &(s->copy_numbers[0]);
591 num_lattice_sites = s->copy_numbers.size();
592 //*positions = s->grid->get_raw_positions().data();
593 *positions = &(s->grid->get_raw_positions()[0]);
594 }
595 return num_lattice_sites;
596 }
597
598