1 /* Steven Andrews, started 10/22/2001.
2  This is a library of functions for the Smoldyn program.
3  See documentation called SmoldynManual.pdf and SmoldynCodeDoc.pdf, and the Smoldyn
4  website, which is at www.smoldyn.org.
5  Copyright 2003-2016 by Steven Andrews.  This work is distributed under the terms
6  of the Gnu Lesser General Public License (LGPL). */
7 
8 #include <stdio.h>
9 #include <math.h>
10 #include <string.h>
11 #include "Geometry.h"
12 #include "math2.h"
13 #include "random2.h"
14 #include "Rn.h"
15 #include "Sphere.h"
16 #include "string2.h"
17 
18 #include "smoldyn.h"
19 #include "smoldynfuncs.h"
20 
21 
22 /******************************************************************************/
23 /********************************** Filaments *********************************/
24 /******************************************************************************/
25 
26 
27 /******************************************************************************/
28 /****************************** Local declarations ****************************/
29 /******************************************************************************/
30 
31 // enumerated types
32 
33 
34 // low level utilities
35 double filStretchEnergy(filamentptr fil,int seg1,int seg2);
36 double filBendEnergy(filamentptr fil,int seg1,int seg2);
37 
38 // memory management
39 void beadfree(beadptr bead);
40 void segmentfree(segmentptr segment);
41 void filamenttypefree(filamenttypeptr filtype);
42 filamentptr filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxmonomer);
43 void filfree(filamentptr fil);
44 
45 
46 // data structure output
47 
48 // structure set up
49 int filtypeSetColor(filamenttypeptr filtype,const double *rgba);
50 int filtypeSetEdgePts(filamenttypeptr filtype,double value);
51 int filtypeSetStipple(filamenttypeptr filtype,int factor,int pattern);
52 int filtypeSetDrawmode(filamenttypeptr filtype,enum DrawMode dm);
53 int filtypeSetShiny(filamenttypeptr filtype,double shiny);
54 int filtypeSetBiology(filamenttypeptr filtype,enum FilamentBiology fb);
55 
56 int filupdateparams(simptr sim);
57 int filupdatelists(simptr sim);
58 
59 
60 
61 /******************************************************************************/
62 /********************************* enumerated types ***************************/
63 /******************************************************************************/
64 
65 
66 /* filFD2string */
filFD2string(enum FilamentDynamics fd,char * string)67 char *filFD2string(enum FilamentDynamics fd,char *string) {
68 	if(fd==FDnone) strcpy(string,"none");
69 	else if(fd==FDrouse) strcpy(string,"rouse");
70 	else if(fd==FDalberts) strcpy(string,"alberts");
71 	else if(fd==FDnedelec) strcpy(string,"nedelec");
72 	else strcpy(string,"none");
73 	return string; }
74 
75 
76 /* filstring2FD */
filstring2FD(const char * string)77 enum FilamentDynamics filstring2FD(const char *string) {
78 	enum FilamentDynamics ans;
79 
80 	if(strbegin(string,"none",0)) ans=FDnone;
81 	else if(strbegin(string,"rouse",0)) ans=FDrouse;
82 	else if(strbegin(string,"alberts",0)) ans=FDalberts;
83 	else if(strbegin(string,"nedelec",0)) ans=FDnedelec;
84 	else ans=FDnone;
85 	return ans; }
86 
87 
88 /* filFB2string */
filFB2string(enum FilamentBiology fb,char * string)89 char *filFB2string(enum FilamentBiology fb,char *string) {
90 	if(fb==FBactin) strcpy(string,"actin");
91 	else if(fb==FBmicrotubule) strcpy(string,"microtubule");
92 	else if(fb==FBintermediate) strcpy(string,"intermediate");
93 	else if(fb==FBdsDNA) strcpy(string,"dsDNA");
94 	else if(fb==FBssDNA) strcpy(string,"ssDNA");
95 	else strcpy(string,"other");
96 	return string; }
97 
98 
99 /* filstring2FB */
filstring2FB(const char * string)100 enum FilamentBiology filstring2FB(const char *string) {
101 	enum FilamentBiology ans;
102 
103 	if(strbegin(string,"actin",0)) ans=FBactin;
104 	else if(strbegin(string,"microtubule",0)) ans=FBmicrotubule;
105 	else if(strbegin(string,"intermediate",0)) ans=FBintermediate;
106 	else if(strbegin(string,"dsDNA",0)) ans=FBdsDNA;
107 	else if(strbegin(string,"ssDNA",0)) ans=FBssDNA;
108 	else if(strbegin(string,"other",0)) ans=FBother;
109 	else ans=FBnone;
110 	return ans; }
111 
112 
113 
114 /******************************************************************************/
115 /****************************** low level utilities ***************************/
116 /******************************************************************************/
117 
118 
119 /* filRandomLength */
filRandomLength(const filamenttypeptr filtype,double thickness,double sigmamult)120 double filRandomLength(const filamenttypeptr filtype,double thickness,double sigmamult) {
121 	double lsigma;
122 	double length;
123 
124 	if(filtype->klen<=0) return filtype->stdlen;
125 	lsigma=sigmamult*sqrt(filtype->kT/(thickness*filtype->klen));
126 	length=0;
127 	while(length<=0)
128 		length=filtype->stdlen+lsigma*gaussrandD();
129 	return length; }
130 
131 
132 /* filRandomAngle */
filRandomAngle(const filamenttypeptr filtype,double thickness,double * angle,double sigmamult)133 double *filRandomAngle(const filamenttypeptr filtype,double thickness,double *angle,double sigmamult) {
134 	double sigma;
135 	int d;
136 
137 	for(d=0;d<3;d++){
138 		if(filtype->kypr[d]>0)
139 			sigma=sigmamult*sqrt(filtype->kT/(thickness*filtype->kypr[d]));
140 		else if(filtype->kypr[d]==0)
141 			sigma=unirandOCD(-PI,PI);
142 		else
143 			sigma=0;
144 		angle[d]=filtype->stdypr[d]+(sigma>0?sigma*gaussrandD():0); }
145 	return angle; }
146 
147 
148 /******************************************************************************/
149 /*************************** computations on filaments ************************/
150 /******************************************************************************/
151 
152 
153 // filStretchEnergy
filStretchEnergy(const filamentptr fil,int seg1,int seg2)154 double filStretchEnergy(const filamentptr fil,int seg1,int seg2) {
155 	double thk,klen,len,stdlen,energy,*xyz,*xyzplus;
156 	filamenttypeptr filtype;
157 	int bs;
158 
159 	filtype=fil->filtype;
160 	stdlen=filtype->stdlen;
161 	klen=filtype->klen;
162 	if(klen<=0) return 0;
163 
164 	energy=0;
165 	if(seg1==-1) seg1=fil->frontbs;
166 	if(seg2==-1) seg2=fil->frontbs+fil->nbs;
167 	if(filtype->isbead) {
168 		for(bs=seg1;bs<seg2-1;bs++) {
169 			xyz=fil->beads[bs]->xyz;
170 			xyzplus=fil->beads[bs+1]->xyz;
171 			len=sqrt((xyzplus[0]-xyz[0])*(xyzplus[0]-xyz[0])+(xyzplus[1]-xyz[1])*(xyzplus[1]-xyz[1])+(xyzplus[2]-xyz[2])*(xyzplus[2]-xyz[2]));
172 			energy+=0.5*klen*(len-stdlen)*(len-stdlen); }}
173 	else {
174 		for(bs=seg1;bs<seg2;bs++){
175 			thk=fil->segments[bs]->thk;
176 			len=fil->segments[bs]->len;
177 			energy+=0.5*thk*klen*(len-stdlen)*(len-stdlen); }}
178 	return energy; }
179 
180 
181 // filBendEnergy
filBendEnergy(const filamentptr fil,int seg1,int seg2)182 double filBendEnergy(const filamentptr fil,int seg1,int seg2) {
183 	double *kypr,*stdypr,*ypr,thk,thkminus,energy;
184 	filamenttypeptr filtype;
185 	segmentptr segptr,segptrminus;
186 	double tk;
187 	int seg;
188 
189 	filtype=fil->filtype;
190 	if(filtype->isbead) return 0;
191 
192 	kypr=filtype->kypr;
193 	stdypr=filtype->stdypr;
194 
195 	energy=0;
196 	if(seg1==-1) seg1=fil->frontbs;
197 	if(seg2==-1) seg2=fil->frontbs+fil->nbs;
198 	for(seg=seg1+1;seg<seg2;seg++) {
199 		segptr = fil->segments[seg];
200 		segptrminus = fil->segments[seg-1];
201 		ypr = segptr->ypr;
202 		thk = segptr->thk;
203 		thkminus = segptrminus->thk;
204 		tk=0.5*(thkminus+thk);
205 		if(kypr[0]>0) energy+=0.5*tk*kypr[0]*(ypr[0]-stdypr[0])*(ypr[0]-stdypr[0]);
206 		if(kypr[1]>0) energy+=0.5*tk*kypr[1]*(ypr[1]-stdypr[1])*(ypr[1]-stdypr[1]);
207 		if(kypr[2]>0) energy+=0.5*tk*kypr[2]*(ypr[2]-stdypr[2])*(ypr[2]-stdypr[2]); }
208 	return energy; }
209 
210 
211 /******************************************************************************/
212 /******************************* memory management ****************************/
213 /******************************************************************************/
214 
215 /* beadalloc */
beadalloc()216 beadptr beadalloc() {
217 	beadptr bead;
218 
219 	CHECKMEM(bead=(beadptr) malloc(sizeof(struct beadstruct)));
220 	bead->xyz[0]=bead->xyz[1]=bead->xyz[2]=0;
221 	bead->xyzold[0]=bead->xyzold[1]=bead->xyzold[2]=0;
222 	return bead;
223 
224  failure:
225 	return NULL; }
226 
227 
228 /* beadfree */
beadfree(beadptr bead)229 void beadfree(beadptr bead) {
230 	if(!bead) return;
231 	free(bead);
232 	return; }
233 
234 
235 /* segmentalloc */
segmentalloc()236 segmentptr segmentalloc() {
237 	segmentptr segment;
238 
239 	CHECKMEM(segment=(segmentptr) malloc(sizeof(struct segmentstruct)));
240 	segment->xyzfront[0]=segment->xyzfront[1]=segment->xyzfront[2]=0;
241 	segment->xyzback[0]=segment->xyzback[1]=segment->xyzback[2]=0;
242 	segment->len=0;
243 	segment->ypr[0]=segment->ypr[1]=segment->ypr[2]=0;
244 	Sph_One2Dcm(segment->dcm);
245 	Sph_One2Dcm(segment->adcm);
246 	segment->thk=1;
247 	return segment;
248 
249  failure:
250 	return NULL; }
251 
252 
253 /* segmentfree */
segmentfree(segmentptr segment)254 void segmentfree(segmentptr segment) {
255 	if(!segment) return;
256 	free(segment);
257 	return; }
258 
259 
260 /* filalloc */
filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxmonomer)261 filamentptr filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxmonomer) {
262 	int bs,br,mn;
263 	beadptr *newbeads;
264 	segmentptr *newsegments;
265 	int *newbranchspots;
266 	filamentptr *newbranches;
267 	char *newmonomers;
268 
269 	if(!fil) {
270 		CHECKMEM(fil=(filamentptr) malloc(sizeof(struct filamentstruct)));
271 		fil->filtype=NULL;
272 		fil->filname=NULL;
273 		fil->maxbs=0;
274 		fil->nbs=0;
275 		fil->frontbs=0;
276 		fil->beads=NULL;
277 		fil->segments=NULL;
278 		fil->frontend=NULL;
279 		fil->backend=NULL;
280 		fil->maxbranch=0;
281 		fil->nbranch=0;
282 		fil->branchspots=NULL;
283 		fil->branches=NULL;
284 		fil->maxmonomer=0;
285 		fil->nmonomer=0;
286 		fil->frontmonomer=0;
287 		fil->monomers=NULL; }
288 
289 	if(maxbs>fil->maxbs) {
290 		if(isbead) {
291 			CHECKMEM(newbeads=(beadptr*) calloc(maxbs,sizeof(struct beadstruct)));
292 			for(bs=0;bs<maxbs;bs++)
293 				newbeads[bs]=NULL;
294 			for(bs=0;bs<fil->maxbs;bs++)
295 				newbeads[bs]=fil->beads[bs];
296 			for(;bs<maxbs;bs++)
297 				CHECKMEM(newbeads[bs]=beadalloc());
298 			free(fil->beads);
299 			fil->beads=newbeads; }
300 		else {
301 			CHECKMEM(newsegments=(segmentptr*) calloc(maxbs,sizeof(struct segmentstruct)));
302 			for(bs=0;bs<maxbs;bs++)
303 				newsegments[bs]=NULL;
304 			for(bs=0;bs<fil->maxbs;bs++)
305 				newsegments[bs]=fil->segments[bs];
306 			for(;bs<maxbs;bs++)
307 				CHECKMEM(newsegments[bs]=segmentalloc());
308 			free(fil->segments);
309 			fil->segments=newsegments; }
310 		fil->maxbs=maxbs; }
311 
312 	if(maxbranch>fil->maxbranch) {
313 		CHECKMEM(newbranchspots=(int*) calloc(maxbranch,sizeof(int)));
314 		CHECKMEM(newbranches=(filamentptr*) calloc(maxbranch,sizeof(filamentptr)));
315 		for(br=0;br<fil->maxbranch;br++) {
316 			newbranchspots[br]=fil->branchspots[br];
317 			newbranches[br]=fil->branches[br]; }
318 		for(;br<maxbranch;br++) {
319 			newbranchspots[br]=0;
320 			newbranches[br]=NULL; }
321 		free(fil->branches);
322 		fil->branches=newbranches;
323 		fil->maxbranch=maxbranch;	}
324 
325 	if(maxmonomer>fil->maxmonomer) {
326 		CHECKMEM(newmonomers=(char*) calloc(maxmonomer,sizeof(char)));
327 		for(mn=0;mn<fil->maxmonomer;mn++)
328 			newmonomers[mn]=fil->monomers[mn];
329 		for(;mn<maxmonomer;mn++)
330 			newmonomers[mn]='\0';
331 		free(fil->monomers);
332 		fil->monomers=newmonomers;
333 		fil->maxmonomer=maxmonomer; }
334 
335 	return fil;
336 
337 failure:		// no memory is freed, which causes leaks, but otherwise risk segfault
338 	return NULL; }
339 
340 
341 /* filfree */
filfree(filamentptr fil)342 void filfree(filamentptr fil) {
343 	int bs;
344 
345 	if(!fil) return;
346 	if(fil->beads) {
347 		for(bs=0;bs<fil->maxbs;bs++)
348 			beadfree(fil->beads[bs]);
349 		free(fil->beads); }
350 	if(fil->segments) {
351 		for(bs=0;bs<fil->maxbs;bs++)
352 			segmentfree(fil->segments[bs]);
353 		free(fil->segments); }
354 
355 	free(fil->branchspots);
356 	free(fil->branches);
357 	free(fil->monomers);
358 	free(fil);
359 	return; }
360 
361 
362 /* filamemttypealloc */
filamenttypealloc(filamenttypeptr filtype,int maxfil,int maxface)363 filamenttypeptr filamenttypealloc(filamenttypeptr filtype,int maxfil,int maxface) {
364 	int f,fc;
365 	filamentptr *newfillist;
366 	char **newfacename,**newfilnames;
367 
368 	if(!filtype) {
369 		CHECKMEM(filtype=(filamenttypeptr) malloc(sizeof(struct filamenttypestruct)));
370 		filtype->filss=NULL;
371 		filtype->ftname=NULL;
372 		filtype->dynamics=FDnone;
373 		filtype->isbead=0;
374 		filtype->biology=FBother;
375 		filtype->bundlevalue=1;
376 
377 		filtype->color[0]=filtype->color[1]=filtype->color[2]=0;
378 		filtype->color[3]=1;
379 		filtype->edgepts=1;
380 		filtype->edgestipple[0]=1;
381 		filtype->edgestipple[1]=0xFFFF;
382 		filtype->drawmode=DMedge;
383 		filtype->shiny=0;
384 
385 		filtype->stdlen=1;
386 		filtype->stdypr[0]=filtype->stdypr[1]=filtype->stdypr[2]=0;
387 		filtype->klen=1;
388 		filtype->kypr[0]=filtype->kypr[1]=filtype->kypr[2]=1;
389 		filtype->kT=0;
390 		filtype->treadrate=0;
391 		filtype->viscosity=1;
392 		filtype->filradius=1;
393 
394 		filtype->maxface=0;
395 		filtype->nface=0;
396 		filtype->facename=NULL;
397 		filtype->facetwist=0;
398 
399 		filtype->maxfil=0;
400 		filtype->nfil=0;
401 		filtype->fillist=NULL;
402 		filtype->filnames=NULL; }
403 
404 	if(maxfil>filtype->maxfil) {
405 		CHECKMEM(newfillist=(filamentptr*) calloc(maxfil,sizeof(filamentptr)));
406 		CHECKMEM(newfilnames=(char**) calloc(maxfil,sizeof(char*)));
407 		for(f=0;f<maxfil;f++) {
408 			newfillist[f]=NULL;
409 			newfilnames[f]=NULL; }
410 		for(f=0;f<filtype->maxfil;f++) {
411 			newfillist[f]=filtype->fillist[f];
412 			newfilnames[f]=filtype->filnames[f]; }
413 		for(;f<maxfil;f++) {
414 			CHECKMEM(newfillist[f]=filalloc(NULL,0,0,0,0));
415 			CHECKMEM(newfilnames[f]=EmptyString());
416 			newfillist[f]->filtype=filtype;
417 			newfillist[f]->filname=newfilnames[f]; }
418 		free(filtype->fillist);
419 		free(filtype->filnames);
420 		filtype->fillist=newfillist;
421 		filtype->filnames=newfilnames;
422 		filtype->maxfil=maxfil; }
423 
424 	if(maxface>filtype->maxface) {
425 		CHECKMEM(newfacename=(char**) calloc(maxface,sizeof(char *)));
426 		for(fc=0;fc<maxface;fc++)
427 			newfacename[fc]=NULL;
428 		for(fc=0;fc<filtype->maxface;fc++)
429 			newfacename[fc]=filtype->facename[fc];
430 		for(;fc<maxface;fc++)
431 			CHECKMEM(newfacename[fc]=EmptyString());
432 		filtype->facename=newfacename;
433 		filtype->maxface=maxface; }
434 
435 	return filtype;
436 
437  failure:
438 	return NULL; }
439 
440 
441 /* filamenttypefree */
filamenttypefree(filamenttypeptr filtype)442 void filamenttypefree(filamenttypeptr filtype) {
443 	if(!filtype) return;
444 	int f,fc;
445 
446 	if(filtype->fillist && filtype->filnames) {
447 		for(f=0;f<filtype->maxfil;f++) {
448 			filfree(filtype->fillist[f]);
449 			free(filtype->filnames[f]); }
450 		free(filtype->fillist);
451 		free(filtype->filnames); }
452 
453 	if(filtype->facename) {
454 		for(fc=0;fc<filtype->maxface;fc++)
455 			free(filtype->facename[fc]);
456 		free(filtype->facename); }
457 
458 	free(filtype);
459 	return; }
460 
461 
462 /* filssalloc */
filssalloc(filamentssptr filss,int maxtype)463 filamentssptr filssalloc(filamentssptr filss,int maxtype) {
464 	char **newftnames;
465 	filamenttypeptr *newfiltypes;
466 	int ft;
467 
468 	if(!filss) {
469 		CHECKMEM(filss=(filamentssptr) malloc(sizeof(struct filamentsuperstruct)));
470 		filss->condition=SCinit;
471 		filss->sim=NULL;
472 		filss->maxtype=0;
473 		filss->ntype=0;
474 		filss->ftnames=NULL;
475 		filss->filtypes=NULL; }
476 
477 	if(maxtype>=filss->maxtype) {
478 		CHECKMEM(newfiltypes=(filamenttypeptr*) calloc(maxtype,sizeof(filamenttypeptr)));
479 		CHECKMEM(newftnames=(char**) calloc(maxtype,sizeof(char*)));
480 		for(ft=0;ft<maxtype;ft++) {
481 			newfiltypes[ft]=NULL;
482 			newftnames[ft]=NULL; }
483 		for(ft=0;ft<filss->maxtype;ft++) {
484 			newfiltypes[ft]=filss->filtypes[ft];
485 			newftnames[ft]=filss->ftnames[ft]; }
486 		for(;ft<maxtype;ft++) {
487 			CHECKMEM(newfiltypes[ft]=filamenttypealloc(NULL,0,0));
488 			CHECKMEM(newftnames[ft]=EmptyString());
489 			newfiltypes[ft]->filss=filss;
490 			newfiltypes[ft]->ftname=newftnames[ft]; }
491 
492 		free(filss->ftnames);
493 		filss->ftnames=newftnames;
494 		free(filss->filtypes);
495 		filss->filtypes=newfiltypes;
496 		filss->maxtype=maxtype; }
497 
498 	return filss;
499 
500 failure:
501 	return NULL; }
502 
503 
504 /* filssfree */
filssfree(filamentssptr filss)505 void filssfree(filamentssptr filss) {
506 	int ft;
507 
508 	if(!filss) return;
509 
510 	if(filss->filtypes) {
511 		for(ft=0;ft<filss->maxtype;ft++)
512 			filamenttypefree(filss->filtypes[ft]);
513 		free(filss->filtypes); }
514 
515 	if(filss->ftnames) {
516 		for(ft=0;ft<filss->maxtype;ft++)
517 			free(filss->ftnames[ft]);
518 		free(filss->ftnames); }
519 
520 	free(filss);
521 	return; }
522 
523 
524 /******************************************************************************/
525 /***************************** data structure output **************************/
526 /******************************************************************************/
527 
528 
529 /* filoutput */
filoutput(const filamentptr fil)530 void filoutput(const filamentptr fil) {
531 	int bs,br,mn,isbead,dim;
532 	simptr sim;
533 	beadptr bead;
534 	segmentptr segment;
535 
536 	if(!fil) {
537 		simLog(NULL,2," NULL filament\n");
538 		return; }
539 
540 	isbead=fil->filtype?fil->filtype->isbead:0;
541 	if(fil->filtype && fil->filtype->filss) {
542 		sim=fil->filtype->filss->sim;
543 		dim=sim->dim; }
544 	else {
545 		sim=NULL;
546 		dim=3; }
547 
548 	simLog(sim,2," Filament: %s\n",fil->filname);
549 	simLog(sim,1,"  allocated beads or segments: %i\n",fil->maxbs);
550 	simLog(sim,2,"  number of %s: %i\n",isbead?"beads":"segments",fil->nbs);
551 	simLog(sim,1,"  front index: %i\n",fil->frontbs);
552 
553 	if(isbead)
554 		simLog(sim,2,"  bead, position\n");
555 	else
556 		simLog(sim,2,"  segment, thickness, length, front position, relative angle\n");
557 	for(bs=0;bs<fil->nbs;bs++) {
558 		if(isbead) {
559 			bead=fil->beads[bs+fil->frontbs];
560 			if(dim==2)
561 				simLog(sim,bs>5?1:2,"   %i pos.=(%1.3g %1.3g)\n",bs,bead->xyz[0],bead->xyz[1]);
562 			else
563 				simLog(sim,bs>5?1:2,"   %i x=(%1.3g %1.3g %1.3g)\n",bs,bead->xyz[0],bead->xyz[1],bead->xyz[2]); }
564 		else {
565 			segment=fil->segments[bs+fil->frontbs];
566 			if(dim==2)
567 				simLog(sim,bs>5?1:2,"   %i thick=%1.3g, length=%1.3g, front pos.=(%1.3g %1.3g), rel. angle=%1.3g\n",bs,segment->thk,segment->len,segment->xyzfront[0],segment->xyzfront[1],segment->ypr[0]);
568 			else
569 				simLog(sim,bs>5?1:2,"   %i thick=%1.3g, length=%1.3g, front pos.=(%1.3g %1.3g %1.3g), rel. angle=(%1.3g %1.3g %1.3g)\n",bs,segment->thk,segment->len,segment->xyzfront[0],segment->xyzfront[1],segment->xyzfront[2],segment->ypr[0],segment->ypr[1],segment->ypr[2]); }}
570 
571 	if(fil->frontend)
572 		simLog(sim,2,"  front branched from: %s\n",fil->frontend->filname);
573 	if(fil->backend)
574 		simLog(sim,2,"  back branched from: %s\n",fil->backend->filname);
575 	if(fil->nbranch>0) {
576 		simLog(sim,2,"  number of branches: %i\n",fil->nbranch);
577 		for(br=0;br<fil->nbranch;br++)
578 			simLog(sim,2,"   %s at %i\n",fil->branches[br]->filname,fil->branchspots[br]); }
579 
580 	if(fil->nmonomer) {
581 		simLog(sim,1,"  monomer codes: %i of %i allocated,",fil->nmonomer,fil->maxmonomer);
582 		simLog(sim,1," starting index: %i\n",fil->frontmonomer);
583 		simLog(sim,2,"  monomer code: ");
584 		for(mn=0;mn<fil->nmonomer;mn++)
585 			simLog(sim,2,"%c",fil->monomers[mn]);
586 		simLog(sim,2,"\n"); }
587 
588 	if(fil->filtype->klen>0)
589 		simLog(sim,2," stretching energy: %g\n",filStretchEnergy(fil,-1,-1));
590 	if(fil->filtype->kypr[0]>0 || fil->filtype->kypr[1]>0 || fil->filtype->kypr[2]>0)
591 		simLog(sim,2," bending energy: %g\n",filBendEnergy(fil,-1,-1));
592 
593 	return; }
594 
595 
596 /* filtypeoutput */
filtypeoutput(const filamenttypeptr filtype)597 void filtypeoutput(const filamenttypeptr filtype) {
598 	char string[STRCHAR];
599 	int isbead,fc,f,dim;
600 	simptr sim;
601 
602 	if(!filtype) {
603 		simLog(NULL,2," NULL filament type\n");
604 		return; }
605 
606 	if(filtype->filss) {
607 		sim=filtype->filss->sim;
608 		dim=sim->dim; }
609 	else {
610 		sim=NULL;
611 		dim=3; }
612 	isbead=filtype->isbead;
613 
614 	simLog(sim,2," Filament type: %s\n",filtype->ftname);
615 	simLog(sim,2,"  dynamics: %s using %s\n",filFD2string(filtype->dynamics,string),isbead?"beads":"segments");
616 	simLog(sim,2,"  biology: %s\n",filFB2string(filtype->biology,string));
617 	simLog(sim,filtype->bundlevalue!=1?2:1,"  bundle value: %g\n",filtype->bundlevalue);
618 
619 	simLog(sim,2,"  color: %g %g %g %g\n",filtype->color[0],filtype->color[1],filtype->color[2],filtype->color[3]);
620 	simLog(sim,2,"  edge points: %g, polygon mode: %s\n",filtype->edgepts,surfdm2string(filtype->drawmode,string));
621 	if(filtype->edgestipple[1]!=0xFFFF) simLog(sim,2,"  edge stippling: %ui %X\n",filtype->edgestipple[0],filtype->edgestipple[1]);
622 	if(filtype->shiny!=0) simLog(sim,2,"  shininess: %g\n",filtype->shiny);
623 
624 	simLog(sim,2,"  %s length: %g\n",filtype->klen>0?"standard":"fixed",filtype->stdlen);
625 	if(filtype->klen>0) simLog(sim,2,"  length force constant: %g\n",filtype->klen);
626 
627 	if(!isbead) {
628 		if(dim==2) {
629 			simLog(sim,2,"  standard angle: %g\n",filtype->stdypr[0]);
630 			simLog(sim,2,"  bending force constant: %g\n",filtype->kypr[0]); }
631 		else {
632 			simLog(sim,2,"  standard angles: %g, %g, %g\n",filtype->stdypr[0],filtype->stdypr[1],filtype->stdypr[2]);
633 			simLog(sim,2,"  bending force constants: %g, %g, %g\n",filtype->kypr[0],filtype->kypr[1],filtype->kypr[2]); }}
634 
635 	simLog(sim,2,"  kT: %g\n",filtype->kT);
636 	simLog(sim,filtype->treadrate>0?2:1,"  treadmilling rate: %g\n",filtype->treadrate);
637 	simLog(sim,2,"  viscosity: %g\n",filtype->viscosity);
638 	simLog(sim,2,"  %s radius: %g\n",isbead?"bead":"segment",filtype->filradius);
639 
640 	if(filtype->nface>0) {
641 		simLog(sim,2,"  %i faces with twist of %g:",filtype->nface,filtype->facetwist);
642 		for(fc=0;fc<filtype->nface;fc++)
643 			simLog(sim,2," %s",filtype->facename[fc]);
644 		simLog(sim,2,"\n"); }
645 
646 	simLog(sim,2,"  %i filaments:\n",filtype->nfil);
647 	for(f=0;f<filtype->nfil;f++)
648 		filoutput(filtype->fillist[f]);
649 
650 	return; }
651 
652 
653 /* filssoutput */
filssoutput(const simptr sim)654 void filssoutput(const simptr sim) {
655 	char string[STRCHAR];
656 	filamentssptr filss;
657 	int ft;
658 
659 	filss=sim->filss;
660 	if(!filss) return;
661 
662 	simLog(sim,2,"FILAMENT PARAMETERS\n");
663 	simLog(sim,filss->condition<SCok?2:1," condition: %s\n",simsc2string(filss->condition,string));
664 	simLog(sim,1," filament types allocated: %i,",filss->maxtype);
665 	simLog(sim,2," filament types defined: %i\n",filss->ntype);
666 	for(ft=0;ft<filss->ntype;ft++)
667 		filtypeoutput(filss->filtypes[ft]);
668 
669 	simLog(sim,2,"\n");
670 	return; }
671 
672 
673 /* filwrite */
filwrite(const simptr sim,FILE * fptr)674 void filwrite(const simptr sim,FILE *fptr) {
675 	filamentssptr filss;
676 
677 	filss=sim->filss;
678 	if(!filss) return;
679 	fprintf(fptr,"# filament parameters\n");
680 	fprintf(fptr,"# ERROR: code not written yet\n\n");
681 	//??TODO: write filwrite
682 	return; }
683 
684 
685 /* filcheckparams */
filcheckparams(const simptr sim,int * warnptr)686 int filcheckparams(const simptr sim,int *warnptr) {
687 	int error,warn;
688 	filamentssptr filss;
689 	char string[STRCHAR];
690 
691 	error=warn=0;
692 	filss=sim->filss;
693 	if(!filss) {
694 		if(warnptr) *warnptr=warn;
695 		return 0; }
696 
697 	if(filss->condition!=SCok) {
698 		warn++;
699 		simLog(sim,7," WARNING: filament structure %s\n",simsc2string(filss->condition,string)); }
700 
701 	//TODO: write checkparams
702 
703 	if(warnptr) *warnptr=warn;
704 	return error; }
705 
706 
707 /******************************************************************************/
708 /****************************** filament manipulation *************************/
709 /******************************************************************************/
710 
711 
712 /* filArrayShift */
filArrayShift(filamentptr fil,int shift)713 void filArrayShift(filamentptr fil,int shift) {
714 	int i,frontbs,backbs;
715 	int isbead;
716 	beadptr nullbead;
717 	segmentptr nullsegment;
718 
719 	isbead=fil->filtype->isbead;
720 	if(!shift)
721 		shift=(fil->maxbs-fil->nbs)/2-fil->frontbs;
722 
723 	frontbs=fil->frontbs;
724 	backbs=fil->frontbs+fil->nbs;
725 
726 	if(shift>0) {
727 		if(backbs+shift>fil->maxbs) shift=fil->maxbs-backbs;
728 		if(isbead) {
729 			for(i=backbs+shift-1;i>=frontbs+shift;i--) {		// i is destination
730 				nullbead = fil->beads[i];
731 				fil->beads[i] = fil->beads[i-shift];
732 				fil->beads[i-shift] = nullbead; }}
733 		else {
734 			for(i=backbs+shift-1;i>=frontbs+shift;i--) {
735 				nullsegment = fil->segments[i];
736 				fil->segments[i]=fil->segments[i-shift];
737 				fil->segments[i-shift] = nullsegment; }}
738 		fil->frontbs+=shift; }
739 
740 	else if(shift<0) {
741 		shift=-shift;								// now shift is positive
742 		if(frontbs-shift<0) shift=frontbs;
743 		if(isbead) {
744 			for(i=frontbs-shift;i<backbs-shift;i++) {				// i is destination
745 				nullbead = fil->beads[i];
746 				fil->beads[i]=fil->beads[i+shift];
747 				fil->beads[i+shift] = nullbead;}}
748 		else {
749 			for(i=frontbs-shift;i<backbs-shift;i++) {
750 				nullsegment = fil->segments[i];
751 				fil->segments[i]=fil->segments[i+shift];
752 				fil->segments[i+shift] = nullsegment; }}
753 		fil->frontbs-=shift; }
754 
755 	return; }
756 
757 
758 /* filAddBead */
filAddBead(filamentptr fil,const double * x,double length,char endchar)759 int filAddBead(filamentptr fil,const double *x,double length,char endchar) {
760 	int seg,dim;
761 	double xyz[3];
762 	filamenttypeptr filtype;
763 	beadptr bead,beadminus,beadplus;
764 
765 	filtype=fil->filtype;
766 	dim=filtype->filss->sim->dim;
767 	if(!filtype->isbead) return -2;	// can't add a bead to a non-bead filament type
768 
769 	if(fil->nbs==fil->maxbs) {
770 		fil=filalloc(fil,1,fil->maxbs*2+1,0,0);
771 		if(!fil) return -1; }		// out of memory
772 
773 	if(endchar=='b') {				// add to back end
774 		if(fil->frontbs+fil->nbs==fil->maxbs)
775 			filArrayShift(fil,0);
776 		seg=fil->frontbs+fil->nbs;
777 		bead=fil->beads[seg];
778 		if(fil->nbs==0) {
779 			bead->xyz[0]=bead->xyzold[0]=x[0];
780 			bead->xyz[1]=bead->xyzold[1]=x[1];
781 			bead->xyz[2]=bead->xyzold[2]=x[2]; }		// xyzold = xyz = x
782 		else {
783 			beadminus=fil->beads[seg-1];
784 			if(dim==2) {
785 				circlerandD(xyz,length);
786 				xyz[2]=0; }
787 			else
788 				sphererandCCD(xyz,length,length);
789 			bead->xyz[0]=bead->xyzold[0]=beadminus->xyz[0]+xyz[0];
790 			bead->xyz[1]=bead->xyzold[1]=beadminus->xyz[1]+xyz[1];
791 			bead->xyz[2]=bead->xyzold[2]=beadminus->xyz[2]+xyz[2]; }
792 		fil->nbs++; }
793 
794 	else {									// add to front end
795 		if(fil->frontbs==0) filArrayShift(fil,0);
796 		if(fil->frontbs==0) filArrayShift(fil,1);								// used if nbs=maxbs-1
797 		seg=fil->frontbs-1;
798 		bead=fil->beads[seg];
799 		if(fil->nbs==0) {
800 			bead->xyz[0]=bead->xyzold[0]=x[0];
801 			bead->xyz[1]=bead->xyzold[1]=x[1];
802 			bead->xyz[2]=bead->xyzold[2]=x[2]; }
803 		else {
804 			beadplus=fil->beads[seg+1];
805 			if(dim==2) {
806 				circlerandD(xyz,length);
807 				xyz[2]=0; }
808 			else
809 				sphererandCCD(xyz,length,length);
810 			bead->xyz[0]=bead->xyzold[0]=beadplus->xyz[0]+xyz[0];
811 			bead->xyz[1]=bead->xyzold[1]=beadplus->xyz[1]+xyz[1];
812 			bead->xyz[2]=bead->xyzold[2]=beadplus->xyz[2]+xyz[2]; }
813 		fil->frontbs--;
814 		fil->nbs++; }
815 
816 	return 0; }
817 
818 
819 /* filAddSegment */
filAddSegment(filamentptr fil,const double * x,double length,const double * angle,double thickness,char endchar)820 int filAddSegment(filamentptr fil,const double *x,double length,const double *angle,double thickness,char endchar) {
821 	int seg,dim;
822     UNUSED(dim);
823 	filamenttypeptr filtype;
824 	segmentptr segment,segmentminus,segmentplus;
825 
826 	filtype=fil->filtype;
827 	dim=filtype->filss->sim->dim;
828 	if(filtype->isbead) return -2;	// can't add a segment to a bead filament type
829 
830 	if(fil->nbs==fil->maxbs) {
831 		fil=filalloc(fil,0,fil->maxbs*2+1,0,0);
832 		if(!fil) return -1; }		// out of memory
833 
834 	if(endchar=='b') {
835 		if(fil->frontbs+fil->nbs==fil->maxbs) filArrayShift(fil,0);
836 		seg=fil->frontbs+fil->nbs;
837 		segment=fil->segments[seg];
838 		segment->len=length;
839 		segment->thk=thickness;
840 		Sph_Xyz2Xyz(angle,segment->ypr);												// ypr = angle
841 		Sph_Xyz2Dcm(angle,segment->dcm);												// A = Dcm(angle)
842 		if(fil->nbs==0) {
843 			segment->xyzfront[0]=x[0];														// x_i = input value
844 			segment->xyzfront[1]=x[1];
845 			segment->xyzfront[2]=x[2];
846 			Sph_Dcm2Dcm(segment->dcm,segment->adcm); }						// B = A
847 		else {
848 			segmentminus=fil->segments[seg-1];
849 			segment->xyzfront[0]=segmentminus->xyzback[0];				// x_i = prior end
850 			segment->xyzfront[1]=segmentminus->xyzback[1];
851 			segment->xyzfront[2]=segmentminus->xyzback[2];
852 			Sph_DcmxDcm(segment->dcm,segmentminus->adcm,segment->adcm);	}	// B_i = A_i . B_{i-1}
853 		Sph_DcmtxUnit(segment->adcm,'x',segment->xyzback,segment->xyzfront,segment->len);		// x_{i+1} = x_i + l_i * BT_i . xhat
854 		fil->nbs++; }
855 
856 	else {
857 		if(fil->frontbs==0) filArrayShift(fil,0);
858 		if(fil->frontbs==0) filArrayShift(fil,1);								// used if nbs=maxbs-1
859 		seg=fil->frontbs;
860 		segment=fil->segments[seg];
861 		segment->len=length;
862 		segment->thk=thickness;
863 		if(fil->nbs==0) {
864 			Sph_Xyz2Dcmt(angle,segment->adcm);									// B_0 = Dcmt(angle)
865 			segment->xyzback[0]=x[0];														// back of segment = input value
866 			segment->xyzback[1]=x[1];
867 			segment->xyzback[2]=x[2]; }
868 		else {
869 			segmentplus=fil->segments[seg+1];
870 			segment->xyzback[0]=segmentplus->xyzfront[0];				// back of segment = next front
871 			segment->xyzback[1]=segmentplus->xyzfront[1];
872 			segment->xyzback[2]=segmentplus->xyzfront[2];
873 			Sph_DcmtxDcm(segmentplus->dcm,segmentplus->adcm,segment->adcm); }	// B_i = AT_{i+1} . B_{i+1}
874 		Sph_Dcm2Dcm(segment->adcm,segment->dcm);									// A_i = B_i
875 		Sph_Dcm2Xyz(segment->dcm,segment->ypr);										// a_0 = Xyz(B_0)
876 		Sph_DcmtxUnit(segment->adcm,'x',segment->xyzfront,segment->xyzback,-segment->len);	// x_i = x_{i+1} - l_i * BT_i . xhat
877 		fil->frontbs--;
878 		fil->nbs++; }
879 	return 0; }
880 
881 
882 /* filAddRandomSegments */
filAddRandomSegments(filamentptr fil,int number,const char * xstr,const char * ystr,const char * zstr,double thickness)883 int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char *ystr,const char *zstr,double thickness) {
884 	int i,dim,er;
885 	double f1,pos[3],angle[3],len;
886 	char **varnames;
887 	double *varvalues;
888 	int nvar;
889 	filamenttypeptr filtype;
890 	simptr sim;
891 
892 	filtype=fil->filtype;
893 	sim=filtype->filss->sim;
894 	dim=sim->dim;
895 	varnames=sim->varnames;
896 	varvalues=sim->varvalues;
897 	nvar=sim->nvar;
898 
899 	if(fil->nbs==0) {								// new filament
900 		systemrandpos(sim,pos);
901 		if(!strcmp(xstr,"u"));
902 		else if(strmathsscanf(xstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[0]=f1;
903 		else return 2;
904 		if(!strcmp(ystr,"u"));
905 		else if(strmathsscanf(ystr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[1]=f1;
906 		else return 2;
907 		if(dim==2) pos[2]=0;
908 		else if(!strcmp(zstr,"u"));
909 		else if(strmathsscanf(zstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[2]=f1;
910 		else return 2; }
911 	else
912 		pos[0]=pos[1]=pos[2]=0;
913 
914 	for(i=0;i<number;i++) {
915 		len=filRandomLength(filtype,thickness,1);
916 		filRandomAngle(filtype,thickness,angle,1);
917 		if(dim==2) angle[1]=angle[2]=0;
918 		er=filAddSegment(fil,pos,len,angle,thickness,'b');
919 		if(er) return er; }
920 
921 	return 0; }
922 
923 
924 /* filAddRandomBeads */
filAddRandomBeads(filamentptr fil,int number,const char * xstr,const char * ystr,const char * zstr)925 int filAddRandomBeads(filamentptr fil,int number,const char *xstr,const char *ystr,const char *zstr) {
926 	int i,dim,er;
927 	double f1,pos[3],len;
928 	char **varnames;
929 	double *varvalues;
930 	int nvar;
931 	filamenttypeptr filtype;
932 	simptr sim;
933 
934 	filtype=fil->filtype;
935 	sim=filtype->filss->sim;
936 	dim=sim->dim;
937 	varnames=sim->varnames;
938 	varvalues=sim->varvalues;
939 	nvar=sim->nvar;
940 
941 	if(fil->nbs==0) {
942 		systemrandpos(sim,pos);
943 		if(!strcmp(xstr,"u"));
944 		else if(strmathsscanf(xstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[0]=f1;
945 		else return 2;
946 		if(!strcmp(ystr,"u"));
947 		else if(strmathsscanf(ystr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[1]=f1;
948 		else return 2;
949 		if(dim==2) pos[2]=0;
950 		else if(!strcmp(zstr,"u"));
951 		else if(strmathsscanf(zstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[2]=f1;
952 		else return 2; }
953 	else
954 		pos[0]=pos[1]=pos[2]=0;
955 
956 	for(i=0;i<number;i++) {
957 		len=filRandomLength(filtype,1,1);
958 		er=filAddBead(fil,pos,len,'b');
959 		if(er) return er; }
960 
961 	return 0; }
962 
963 
964 /* filRemoveSegment */
filRemoveSegment(filamentptr fil,char endchar)965 int filRemoveSegment(filamentptr fil,char endchar) {
966 	int seg;
967 	segmentptr segment;
968 
969 	if(fil->nbs==0) return -1;
970 
971 	if(endchar=='b')
972 		fil->nbs--;
973 	else {
974 		seg=++fil->frontbs;
975 		fil->nbs--;
976 		if(!fil->filtype->isbead) {
977 			segment=fil->segments[seg];
978 			Sph_Dcm2Dcm(segment->adcm,segment->dcm);
979 			Sph_Dcm2Xyz(segment->dcm,segment->ypr); }}
980 	return 0; }
981 
982 
983 /* filTranslate */
filTranslate(filamentptr fil,const double * vect,char func)984 void filTranslate(filamentptr fil,const double *vect,char func) {
985 	int seg;
986 	double shift[3];
987 	segmentptr segment;
988 	beadptr bead;
989 
990 	seg=fil->frontbs;
991 	if(func=='=') {
992 		if(fil->filtype->isbead){
993 			bead=fil->beads[seg];
994 			shift[0]=bead->xyz[0]-vect[0];
995 			shift[1]=bead->xyz[1]-vect[1];
996 			shift[2]=bead->xyz[2]-vect[2]; }
997 		else{
998 			segment=fil->segments[seg];
999 			shift[0]=segment->xyzfront[0]-vect[0];
1000 			shift[1]=segment->xyzfront[1]-vect[1];
1001 			shift[2]=segment->xyzfront[2]-vect[2];}}
1002 	else if(func=='-') {
1003 		shift[0]=-vect[0];
1004 		shift[1]=-vect[1];
1005 		shift[2]=-vect[2]; }
1006 	else {
1007 		shift[0]=vect[0];
1008 		shift[1]=vect[1];
1009 		shift[2]=vect[2]; }
1010 
1011 	if(fil->filtype->isbead){
1012 		for(seg=0;seg<fil->nbs;seg++) {
1013 			bead = fil->beads[seg+fil->frontbs];
1014 			bead->xyz[0]+=shift[0];
1015 			bead->xyz[1]+=shift[1];
1016 			bead->xyz[2]+=shift[2];
1017 			bead->xyzold[0]+=shift[0];
1018 			bead->xyzold[1]+=shift[1];
1019 			bead->xyzold[2]+=shift[2]; }}
1020 	else {
1021 		for(seg=0;seg<fil->nbs;seg++) {
1022 			segment=fil->segments[seg+fil->frontbs];
1023 			segment->xyzfront[0]+=shift[0];
1024 			segment->xyzfront[1]+=shift[1];
1025 			segment->xyzfront[2]+=shift[2];
1026 			segment->xyzback[0]+=shift[0];
1027 			segment->xyzback[1]+=shift[1];
1028 			segment->xyzback[2]+=shift[2]; }}
1029 
1030 	return; }
1031 
1032 
1033 /* filRotateVertex */
filRotateVertex(filamentptr fil,int seg,const double * angle,char endchar,char func)1034 void filRotateVertex(filamentptr fil,int seg,const double *angle,char endchar,char func) {
1035 	(void)fil;
1036 	(void)seg;
1037 	(void)angle;
1038 	(void)endchar;
1039 	(void)func;
1040 	return; }
1041 
1042 
1043 /* filLengthenSegment */
filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func)1044 void filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func) {
1045 	(void)fil;
1046 	(void)seg;
1047 	(void)length;
1048 	(void)endchar;
1049 	(void)func;
1050 	return; }
1051 
1052 
1053 /* filReverseFilament */
filReverseFilament(filamentptr fil)1054 void filReverseFilament(filamentptr fil) {
1055 	(void)fil;
1056 	return; }
1057 
1058 
1059 /* filCopyFilament */
filCopyFilament(filamentptr filto,const filamentptr filfrom)1060 int filCopyFilament(filamentptr filto,const filamentptr filfrom) {
1061 	int isbead,i;
1062 	beadptr beadto,beadfrom;
1063 	segmentptr segmentto,segmentfrom;
1064 
1065 	if(!filto || !filfrom) return 2;
1066 
1067 	if(filfrom->filtype) isbead=filfrom->filtype->isbead;
1068 	else if(filfrom->beads) isbead=1;
1069 	else isbead=0;
1070 
1071 	filto->filtype=filfrom->filtype;
1072 
1073 	filto->nbs=0;
1074 	filto->frontbs=0;
1075 	filto->nbranch=0;
1076 	filto->nmonomer=0;
1077 	filto->frontmonomer=0;
1078 	filto=filalloc(filto,isbead,filfrom->nbs,filfrom->nbranch,filfrom->nmonomer);
1079 	if(!filto) return 1;
1080 
1081 	if(isbead) {
1082 		for(i=0;i<filfrom->nbs;i++) {
1083 			beadto=filto->beads[i];
1084 			beadfrom=filfrom->beads[i+filfrom->frontbs];
1085 			copyVD(beadfrom->xyz,beadto->xyz,3);
1086 			copyVD(beadfrom->xyzold,beadto->xyzold,3); }}
1087 	else {
1088 		for(i=0;i<filfrom->nbs;i++) {
1089 			segmentto=filto->segments[i];
1090 			segmentfrom=filfrom->segments[i+filfrom->frontbs];
1091 			copyVD(segmentfrom->xyzfront,segmentto->xyzfront,3);
1092 			copyVD(segmentfrom->xyzback,segmentto->xyzback,3);
1093 			segmentto->len=segmentfrom->len;
1094 			copyVD(segmentfrom->ypr,segmentto->ypr,3);
1095 			copyVD(segmentfrom->dcm,segmentto->dcm,9);
1096 			copyVD(segmentfrom->adcm,segmentto->adcm,9);
1097 			segmentto->thk=segmentfrom->thk; }}
1098 	filto->nbs=filfrom->nbs;
1099 
1100 	filto->frontend=filfrom->frontend;
1101 	filto->backend=filfrom->backend;
1102 	for(i=0;i<filfrom->nbranch;i++) {
1103 		filto->branchspots[i]=filfrom->branchspots[i];
1104 		filto->branches[i]=filfrom->branches[i]; }
1105 	filto->nbranch=filfrom->nbranch;
1106 
1107 	for(i=0;i<filfrom->nmonomer;i++)
1108 		filto->monomers[i]=filfrom->monomers[i+filfrom->frontmonomer];
1109 	filto->nmonomer=filfrom->nmonomer;
1110 
1111 	return 0; }
1112 
1113 
1114 /******************************************************************************/
1115 /********************************* filament type ******************************/
1116 /******************************************************************************/
1117 
1118 
1119 /* filtypeSetParam */
filtypeSetParam(filamenttypeptr filtype,const char * param,int index,double value)1120 int filtypeSetParam(filamenttypeptr filtype,const char *param,int index,double value) {
1121 	int er;
1122 
1123 	er=0;
1124 	if(!strcmp(param,"stdlen")) {
1125 		if(value<=0) er=2;
1126 		else filtype->stdlen=value; }
1127 
1128 	else if(!strcmp(param,"stdypr")) {
1129 		if(index<0)
1130 			filtype->stdypr[0]=filtype->stdypr[1]=filtype->stdypr[2]=value;
1131 		else filtype->stdypr[index]=value; }
1132 
1133 	else if(!strcmp(param,"klen"))
1134 		filtype->klen=value;
1135 
1136 	else if(!strcmp(param,"kypr")) {
1137 		if(index<0) filtype->kypr[0]=filtype->kypr[1]=filtype->kypr[2]=value;
1138 		else filtype->kypr[index]=value; }
1139 
1140 	else if(!strcmp(param,"kT")) {
1141 		if(value<=0) er=2;
1142 		else filtype->kT=value; }
1143 
1144 	else if(!strcmp(param,"treadrate"))
1145 		filtype->treadrate=value;
1146 
1147 	else if(!strcmp(param,"viscosity")) {
1148 		if(value<=0) er=2;
1149 		else filtype->viscosity=value; }
1150 
1151 	else if(!strcmp(param,"beadradius")) {
1152 		if(value<=0) er=2;
1153 		else filtype->filradius=value; }
1154 
1155 	return er; }
1156 
1157 
1158 /* filtypeSetColor */
filtypeSetColor(filamenttypeptr filtype,const double * rgba)1159 int filtypeSetColor(filamenttypeptr filtype,const double *rgba) {
1160 	int col;
1161 
1162 	for(col=0;col<4;col++)
1163 		if(rgba[col]<0 || rgba[col]>1) return 2;
1164 
1165 	for(col=0;col<4;col++) filtype->color[col]=rgba[col];
1166 	return 0; }
1167 
1168 
1169 /* filtypeSetEdgePts */
filtypeSetEdgePts(filamenttypeptr filtype,double value)1170 int filtypeSetEdgePts(filamenttypeptr filtype,double value) {
1171 	if(value<0) return 2;
1172 	filtype->edgepts=value;
1173 	return 0; }
1174 
1175 
1176 /* filtypeSetStipple */
filtypeSetStipple(filamenttypeptr filtype,int factor,int pattern)1177 int filtypeSetStipple(filamenttypeptr filtype,int factor,int pattern) {
1178 	if(factor>=0) {
1179 		if(factor==0) return 2;
1180 		filtype->edgestipple[0]=(unsigned int) factor; }
1181 	if(pattern>=0) {
1182 		if(pattern>0xFFFF) return 2;
1183 		filtype->edgestipple[1]=(unsigned int) pattern; }
1184 	return 0; }
1185 
1186 
1187 /* filtypeSetDrawmode */
filtypeSetDrawmode(filamenttypeptr filtype,enum DrawMode dm)1188 int filtypeSetDrawmode(filamenttypeptr filtype,enum DrawMode dm) {
1189 	if(dm==DMnone) return 2;
1190 	filtype->drawmode=dm;
1191 	return 0; }
1192 
1193 
1194 /* filtypeSetShiny */
filtypeSetShiny(filamenttypeptr filtype,double shiny)1195 int filtypeSetShiny(filamenttypeptr filtype,double shiny) {
1196 	if(!filtype) return 1;
1197 	if(shiny<0 || shiny>128) return 2;
1198 	filtype->shiny=shiny;
1199 	return 0; }
1200 
1201 
1202 /* filtypeSetDynamics */
filtypeSetDynamics(filamenttypeptr filtype,enum FilamentDynamics fd)1203 int filtypeSetDynamics(filamenttypeptr filtype,enum FilamentDynamics fd) {
1204 	filtype->dynamics=fd;			//?? check for pre-existing dynamics and convert if needed
1205 	if(fd==FDrigidbeads) filtype->isbead=1;
1206 	else if(fd==FDrigidsegments) filtype->isbead=0;
1207 	else if(fd==FDrouse) filtype->isbead=1;
1208 	else if(fd==FDalberts) filtype->isbead=0;
1209 	else if(fd==FDnedelec) filtype->isbead=0;
1210 	return 0; }
1211 
1212 
1213 /* filtypeSetBiology */
filtypeSetBiology(filamenttypeptr filtype,enum FilamentBiology fb)1214 int filtypeSetBiology(filamenttypeptr filtype,enum FilamentBiology fb) {
1215 	filtype->biology=fb;	//?? set parameters to match
1216 	return 0; }
1217 
1218 
1219 /******************************************************************************/
1220 /**************************** filament superstructure *************************/
1221 /******************************************************************************/
1222 
1223 
1224 /* filsetcondition */
filsetcondition(filamentssptr filss,enum StructCond cond,int upgrade)1225 void filsetcondition(filamentssptr filss,enum StructCond cond,int upgrade) {
1226 	if(!filss) return;
1227 	if(upgrade==0 && filss->condition>cond) filss->condition=cond;
1228 	else if(upgrade==1 && filss->condition<cond) filss->condition=cond;
1229 	else if(upgrade==2) filss->condition=cond;
1230 	if(filss->sim && filss->condition<filss->sim->condition) {
1231 		cond=filss->condition;
1232 		simsetcondition(filss->sim,cond==SCinit?SClists:cond,0); }
1233 	return; }
1234 
1235 
1236 /* filenablefilaments */
filenablefilaments(simptr sim)1237 int filenablefilaments(simptr sim) {
1238 	filamentssptr filss;
1239 
1240 	if(sim->filss) return 0;
1241 
1242 	filss=filssalloc(sim->filss,1);
1243 	if(!filss) return 1;
1244 	sim->filss=filss;
1245 	filss->sim=sim;
1246 	filsetcondition(sim->filss,SClists,0);
1247 	return 0; }
1248 
1249 
1250 /* filAddFilament */
filAddFilament(filamenttypeptr filtype,filamentptr fil,const char * filname)1251 filamentptr filAddFilament(filamenttypeptr filtype,filamentptr fil,const char *filname) {
1252 	int f;
1253 	filamentptr fil2;
1254 
1255 	if(!filtype && fil)											// no filtype, yes fil
1256 		return fil;
1257 
1258 	else if(!filtype) {											// no filtype, no fil
1259 		fil=filalloc(NULL,0,0,0,0);
1260 		if(!fil) return NULL;
1261 		fil->filname=EmptyString();
1262 		if(!fil->filname) return NULL; }
1263 
1264 	else {																	// yes filtype, yes or no fil
1265 		f=stringfind(filtype->filnames,filtype->nfil,filname);
1266 		if(f<0) {
1267 			if(filtype->nfil==filtype->maxfil) {
1268 				filtype=filamenttypealloc(filtype,filtype->maxfil*2+1,0);
1269 				if(!filtype) return NULL; }
1270 			f=filtype->nfil++;
1271 			strncpy(filtype->filnames[f],filname,STRCHAR-1);
1272 			filtype->filnames[f][STRCHAR-1]='\0';
1273 			fil2=filtype->fillist[f];
1274 			if(fil) {
1275 				filCopyFilament(fil2,fil);
1276 				fil2->filtype=filtype;
1277 				free(fil->filname);
1278 				filfree(fil); }
1279 			fil=fil2;
1280 			filsetcondition(filtype->filss,SClists,0); }
1281 		else
1282 			fil=filtype->fillist[f]; }
1283 
1284 	return fil; }
1285 
1286 
1287 /* filAddFilamentType */
filAddFilamentType(simptr sim,const char * ftname)1288 filamenttypeptr filAddFilamentType(simptr sim,const char *ftname) {
1289 	int ft;
1290 	filamentssptr filss;
1291 	filamenttypeptr filtype;
1292 	filamentssptr fssptr;
1293 
1294 	filss=sim->filss;
1295 
1296 	ft=stringfind(filss->ftnames,filss->ntype,ftname);
1297 	if(ft<0) {
1298 		if(filss->ntype==filss->maxtype) {
1299 			fssptr=filssalloc(filss,filss->maxtype*2+1);
1300 			if(!fssptr) return NULL; }
1301 		ft=filss->ntype++;
1302 		strncpy(filss->ftnames[ft],ftname,STRCHAR-1);
1303 		filss->ftnames[ft][STRCHAR-1]='\0';
1304 		filtype=filss->filtypes[ft];
1305 		filsetcondition(filss,SClists,0); }
1306 	else
1307 		filtype=filss->filtypes[ft];
1308 
1309 	return filtype; }
1310 
1311 
1312 /******************************************************************************/
1313 /********************************** user input ********************************/
1314 /******************************************************************************/
1315 
1316 
1317 /* filtypereadstring */
filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr filtype,const char * word,char * line2)1318 filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr filtype,const char *word,char *line2) {
1319 	char **varnames;
1320 	double *varvalues;
1321 	int nvar;
1322 
1323 	filamentssptr filss;
1324     UNUSED(filss);
1325 	int itct,er,i1,i2;
1326 	char nm[STRCHAR],nm1[STRCHAR];
1327 	double fltv1[9],f1;
1328 	enum DrawMode dm;
1329 	enum FilamentDynamics fd;
1330 	enum FilamentBiology fb;
1331 
1332 	printf("%s\n",word);//?? debug
1333 	filss=sim->filss;
1334 
1335 	varnames=sim->varnames;
1336 	varvalues=sim->varvalues;
1337 	nvar=sim->nvar;
1338 
1339 	if(!strcmp(word,"name")) {								// name
1340 		itct=sscanf(line2,"%s",nm);
1341 		CHECKS(itct==1,"error reading filament type name");
1342 		filtype=filAddFilamentType(sim,nm);
1343 		CHECKS(filtype,"failed to add filament type");
1344 		CHECKS(!strnword(line2,2),"unexpected text following name"); }
1345 
1346 	else if(!strcmp(word,"dynamics")) {				// dynamics
1347 		CHECKS(filtype,"need to enter filament type name before dynamics");
1348 		itct=sscanf(line2,"%s",nm1);
1349 		CHECKS(itct==1,"dynamics options: RigidBeads, RigidSegements, Rouse, Alberts, Nedelec");
1350 		fd=filstring2FD(nm1);
1351 		CHECKS(fd!=FDnone,"dynamics options: RigidBeads, RigidSegements, Rouse, Alberts, Nedelec");
1352 		er=filtypeSetDynamics(filtype,fd); //?? beware of changing between beads and segments
1353 		CHECKS(!er,"BUG: error setting filament dynamics");
1354 		CHECKS(!strnword(line2,2),"unexpected text following dynamics"); }
1355 
1356 	else if(!strcmp(word,"biology")) {				// biology
1357 		itct=sscanf(line2,"%s",nm);
1358 		CHECKS(itct==1,"biology options: actin, microtubule, intermediate, dsDNA, ssDNA, other");
1359 		fb=filstring2FB(nm);
1360 		CHECKS(fb!=FBnone,"biology options: actin, microtubule, intermediate, dsDNA, ssDNA, other");
1361 		er=filtypeSetBiology(filtype,fb);
1362 		CHECKS(!er,"error setting filament biology");
1363 		CHECKS(!strnword(line2,2),"unexpected text following biology"); }
1364 
1365 	else if(!strcmp(word,"color") || !strcmp(word,"colour")) {		// color
1366 		CHECKS(filtype,"need to enter filament type name before color");
1367 		er=graphicsreadcolor(&line2,fltv1);
1368 		CHECKS(er!=3,"color values need to be between 0 and 1");
1369 		CHECKS(er!=4,"color name not recognized");
1370 		CHECKS(er!=6,"alpha values need to be between 0 and 1");
1371 		CHECKS(er==0,"format is either 3 numbers or color name, and then optional alpha value");
1372 		er=filtypeSetColor(filtype,fltv1);
1373 		CHECKS(!er,"BUG: error in filtypeSetColor");
1374 		CHECKS(!line2,"unexpected text following color"); }
1375 
1376 	else if(!strcmp(word,"thickness")) {				// thickness
1377 		CHECKS(filtype,"need to enter filament type name before thickness");
1378 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1379 		CHECKS(itct==1,"thickness value is missing");
1380 		CHECKS(f1>=0,"thickness value needs to be at least 0");
1381 		er=filtypeSetEdgePts(filtype,f1);
1382 		CHECKS(!er,"BUG: error in filtypeSetEdgePts");
1383 		CHECKS(!strnword(line2,2),"unexpected text following thickness"); }
1384 
1385 	else if(!strcmp(word,"stipple")) {					// stipple
1386 		CHECKS(filtype,"need to enter filament type name before stipple");
1387 		itct=strmathsscanf(line2,"%mi %mi",varnames,varvalues,nvar,&i1,&i2);
1388 		CHECKS(itct==2,"stipple format: factor pattern");
1389 		CHECKS(i1>=1,"stipple factor needs to be >=1");
1390 		CHECKS(i2>=0 && i2 <=0xFFFF,"stipple pattern needs to be between 0x00 and 0xFFFF");
1391 		er=filtypeSetStipple(filtype,i1,i2);
1392 		CHECKS(!er,"BUG: error in filtypeSetStipple");
1393 		CHECKS(!strnword(line2,3),"unexpected text following stipple"); }
1394 
1395 	else if(!strcmp(word,"polygon")) {					// polygon
1396 		CHECKS(filtype,"need to enter filament type name before polygon");
1397 		itct=sscanf(line2,"%s",nm1);
1398 		CHECKS(itct==1,"polygon format: drawmode");
1399 		dm=surfstring2dm(nm1);
1400 		CHECKS(dm!=DMnone,"in polygon, drawing mode is not recognized");
1401 		er=filtypeSetDrawmode(filtype,dm);
1402 		CHECKS(!er,"BUG: error in filtypeSetDrawmode");
1403 		CHECKS(!strnword(line2,2),"unexpected text following polygon"); }
1404 
1405 	else if(!strcmp(word,"shininess")) {				// shininess
1406 		CHECKS(filtype,"need to enter filament type name before shininess");
1407 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1408 		CHECKS(itct==1,"shininess format: value");
1409 		CHECKS(f1>=0 && f1<=128,"shininess value needs to be between 0 and 128");
1410 		er=filtypeSetShiny(filtype,f1);
1411 		CHECKS(!er,"BUG: error in filtypeSetShiny");
1412 		CHECKS(!strnword(line2,2),"unexpected text following shininess"); }
1413 
1414 	else if(!strcmp(word,"kT")) {								// kT
1415 		CHECKS(filtype,"need to enter filament type name before kT");
1416 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1417 		CHECKS(itct==1,"kT format: value");
1418 		CHECKS(f1>0,"kT value needs to be greater than 0");
1419 		filtypeSetParam(filtype,"kT",0,f1);
1420 		CHECKS(!strnword(line2,2),"unexpected text following kT"); }
1421 
1422 	else if(!strcmp(word,"treadmill_rate")) {		// treadmill_rate
1423 		CHECKS(filtype,"need to enter filament type name before treadmill_rate");
1424 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1425 		CHECKS(itct==1,"treadmill_rate format: value");
1426 		filtypeSetParam(filtype,"treadrate",0,f1);
1427 		CHECKS(!strnword(line2,2),"unexpected text following treadrate"); }
1428 
1429 	else if(!strcmp(word,"viscosity")) {			// viscosity
1430 		CHECKS(filtype,"need to enter filament type name before viscosity");
1431 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1432 		CHECKS(itct==1,"viscosity format: value");
1433 		CHECKS(f1>0,"viscosity value needs to be greater than 0");
1434 		filtypeSetParam(filtype,"viscosity",0,f1);
1435 		CHECKS(!strnword(line2,2),"unexpected text following viscosity"); }
1436 
1437 	else if(!strcmp(word,"bead_radius")) {			// beadradius
1438 		CHECKS(filtype,"need to enter filament type name before bead radius");
1439 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1440 		CHECKS(itct==1,"bead_radius format: value");
1441 		filtypeSetParam(filtype,"beadradius",0,f1);
1442 		CHECKS(!strnword(line2,2),"unexpected text following bead radius"); }
1443 
1444 	else if(!strcmp(word,"standard_length")) {	// standard_length
1445 		CHECKS(filtype,"need to enter filament type name before standard_length");
1446 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1447 		CHECKS(itct==1,"standard_length format: value");
1448 		CHECKS(f1>0,"standard_length value needs to be greater than 0");
1449 		filtypeSetParam(filtype,"stdlen",0,f1);
1450 		CHECKS(!strnword(line2,2),"unexpected text following standard_length"); }
1451 
1452 	else if(!strcmp(word,"standard_angle")) {	// standard_angle
1453 		CHECKS(filtype,"need to enter filament type name before standard_angle");
1454 		itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2]);
1455 		CHECKS(itct==3,"standard_angle format: yaw pitch roll");
1456 		filtypeSetParam(filtype,"stdypr",0,fltv1[0]);
1457 		filtypeSetParam(filtype,"stdypr",1,fltv1[1]);
1458 		filtypeSetParam(filtype,"stdypr",2,fltv1[2]);
1459 		CHECKS(!strnword(line2,4),"unexpected text following standard_angle"); }
1460 
1461 	else if(!strcmp(word,"force_length")) {		// force_length
1462 		CHECKS(filtype,"need to enter filament type name before force_length");
1463 		itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1);
1464 		CHECKS(itct==1,"force_length format: value");
1465 		filtypeSetParam(filtype,"klen",0,f1);
1466 		CHECKS(!strnword(line2,2),"unexpected text following force_length"); }
1467 
1468 	else if(!strcmp(word,"force_angle")) {		// force_angle
1469 		CHECKS(filtype,"need to enter filament type name before force_angle");
1470 		itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2]);
1471 		CHECKS(itct==3,"force_angle format: yaw pitch roll");
1472 		filtypeSetParam(filtype,"kypr",0,fltv1[0]);
1473 		filtypeSetParam(filtype,"kypr",1,fltv1[1]);
1474 		filtypeSetParam(filtype,"kypr",2,fltv1[2]);
1475 		CHECKS(!strnword(line2,4),"unexpected text following force_angle"); }
1476 
1477 //?? Need to add faces and facetwist
1478 
1479 	else {																				// unknown word
1480 		CHECKS(0,"syntax error within filament_type block: statement not recognized"); }
1481 
1482 	return filtype;
1483 
1484 failure:
1485 	simParseError(sim,pfp);
1486 	return NULL; }
1487 
1488 
1489 /* filreadstring */
filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamenttypeptr filtype,const char * word,char * line2)1490 filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamenttypeptr filtype,const char *word,char *line2) {
1491 	char **varnames;
1492 	double *varvalues;
1493 	int nvar;
1494 
1495 	filamentptr fil2;
1496 	filamentssptr filss;
1497 	int itct,er,i1;
1498 	char nm[STRCHAR],nm1[STRCHAR],endchar,str1[STRCHAR],str2[STRCHAR],str3[STRCHAR],symbol;
1499 	double fltv1[9],length,angle[3],thick;
1500 
1501 	printf("%s\n",word);//?? debug
1502 	filss=sim->filss;
1503 
1504 	varnames=sim->varnames;
1505 	varvalues=sim->varvalues;
1506 	nvar=sim->nvar;
1507 
1508 	if(!strcmp(word,"name")) {								// name
1509 		itct=sscanf(line2,"%s",nm);
1510 		CHECKS(itct==1,"error reading filament name");
1511 		fil=filAddFilament(filtype,fil,nm);
1512 		CHECKS(fil,"failed to add filament");
1513 		CHECKS(!strnword(line2,2),"unexpected text following filament name"); }
1514 
1515 	else if(!strcmp(word,"type")) {						// type
1516 		itct=sscanf(line2,"%s",nm1);
1517 		CHECKS(itct==1,"error reading filament type");
1518 		i1=stringfind(filss->ftnames,filss->ntype,nm1);
1519 		CHECKS(i1>=0,"filament type not recognized");
1520 		filtype=filss->filtypes[i1];
1521 		if(fil) {
1522 			CHECKS(fil->filtype==NULL || fil->filtype==filtype,"filament type was already defined and cannot be changed");
1523 			fil=filAddFilament(filtype,fil,fil->filname);
1524 			CHECKS(fil,"failed to add filament"); }
1525 		CHECKS(!strnword(line2,2),"unexpected text following name"); }
1526 
1527 	else if(!strcmp(word,"first_segment")) {		// first_segment
1528 		CHECKS(fil,"need to enter filament name before first_segment");
1529 		CHECKS(fil->nbs==0,"filament already has segments in it");
1530 		itct=strmathsscanf(line2,"%mlg %mlg %mlg %mlg %mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2],&length,&angle[0],&angle[1],&angle[2]);
1531 		CHECKS(itct==7,"first_segment format: x y z length angle0 angle1 angle2 [thickness]");
1532 		CHECKS(length>0,"length needs to be >0");
1533 		line2=strnword(line2,8);
1534 		thick=1;
1535 		if(line2) {
1536 			itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&thick);
1537 			CHECKS(itct==1,"first_segment format: x y z length angle0 angle1 angle2 [thickness]");
1538 			CHECKS(thick>0,"thickness needs to be >0");
1539 			line2=strnword(line2,2); }
1540 		er=filAddSegment(fil,fltv1,length,angle,thick,'b');
1541 		CHECKS(er==0,"BUG: error in filAddsegment");
1542 		CHECKS(!line2,"unexpected text following first_segment"); }
1543 
1544 	else if(!strcmp(word,"add_segment")) {			// add_segment
1545 		CHECKS(fil,"need to enter filament name before add_segment");
1546 		CHECKS(fil->nbs>0,"use first_segment to enter the first segment");
1547 		itct=strmathsscanf(line2,"%mlg %mlg %mlg %mlg",varnames,varvalues,nvar,&length,&angle[0],&angle[1],&angle[2]);
1548 		CHECKS(itct==4,"add_segment format: length angle0 angle1 angle2 [thickness [end]]");
1549 		CHECKS(length>0,"length needs to be >0");
1550 		line2=strnword(line2,5);
1551 		thick=1;
1552 		endchar='b';
1553 		if(line2) {
1554 			itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&thick);
1555 			CHECKS(itct==1,"add_segment format: length angle0 angle1 angle2 [thickness [end]]");
1556 			CHECKS(thick>0,"thickness needs to be >0");
1557 			line2=strnword(line2,2);
1558 			if(line2) {
1559 				itct=sscanf(line2,"%s",nm1);
1560 				CHECKS(itct==1,"add_segment format: length angle0 angle1 angle2 [thickness [end]]");
1561 				if(nm1[0]=='B' || nm1[0]=='b') endchar='b';
1562 				else if(nm1[0]=='F' || nm1[0]=='f') endchar='f';
1563 				else CHECKS(0,"end needs to be 'back' or 'front'");
1564 				line2=strnword(line2,2); }}
1565 		er=filAddSegment(fil,NULL,length,angle,thick,endchar);
1566 		CHECKS(er==0,"BUG: error in filAddsegment");
1567 		CHECKS(!line2,"unexpected text following add_segment"); }
1568 
1569 	else if(!strcmp(word,"remove_segment")) {			// remove_segment
1570 		CHECKS(fil,"need to enter filament name before remove_segment");
1571 		CHECKS(fil->nbs>0,"filament has no segments to remove");
1572 		itct=sscanf(line2,"%s",nm1);
1573 		CHECKS(itct==1,"remove_segment format: end");
1574 		endchar='b';
1575 		if(nm1[0]=='B' || nm1[0]=='b') endchar='b';
1576 		else if(nm1[0]=='F' || nm1[0]=='f') endchar='f';
1577 		else CHECKS(0,"end needs to be 'back' or 'front'");
1578 		er=filRemoveSegment(fil,endchar);
1579 		CHECKS(!er,"BUG: error in filRemovesegment");
1580 		CHECKS(!strnword(line2,2),"unexpected text following remove_segment"); }
1581 
1582 	else if(!strcmp(word,"random_segments")) {		// random_segments
1583 		CHECKS(fil,"need to enter filament name before random_segments");
1584 		CHECKS(fil->filtype,"need to enter filament type before random_segments");
1585 		itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1586 		CHECKS(itct==1,"random_segments format: number [x y z] [thickness]");
1587 		CHECKS(i1>0,"number needs to be >0");
1588 		line2=strnword(line2,2);
1589 		if(fil->nbs==0) {
1590 			CHECKS(line2,"missing x y z position information");
1591 			itct=sscanf(line2,"%s %s %s",str1,str2,str3);
1592 			CHECKS(itct==3,"random_segments format: number [x y z] [thickness]");
1593 			line2=strnword(line2,4); }
1594 		else {
1595 			sprintf(str1,"%i",0);
1596 			sprintf(str2,"%i",0);
1597 			sprintf(str3,"%i",0); }
1598 		thick=1;
1599 		if(line2) {
1600 			itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&thick);
1601 			CHECKS(itct==1,"random_segments format: number [x y z] [thickness]");
1602 			CHECKS(thick>0,"thickness needs to be >0");
1603 			line2=strnword(line2,2); }
1604 		er=filAddRandomSegments(fil,i1,str1,str2,str3,thick);
1605 		CHECKS(er!=2,"random_segments positions need to be 'u' or value");
1606 		CHECKS(er==0,"BUG: error in filAddRandomsegments");
1607 		CHECKS(!line2,"unexpected text following random_segments"); }
1608 
1609 	else if(!strcmp(word,"translate")) {					// translate
1610 		CHECKS(fil,"need to enter filament name before translate");
1611 		itct=strmathsscanf(line2,"%c %mlg %mlg %mlg",varnames,varvalues,nvar,&symbol,&fltv1[0],&fltv1[1],&fltv1[2]);
1612 		CHECKS(itct==4,"translate format: symbol x y z");
1613 		CHECKS(symbol=='=' || symbol=='+' || symbol=='-',"symbol needs to be '=', '+', or '-'");
1614 		filTranslate(fil,fltv1,symbol);
1615 		CHECKS(!strnword(line2,5),"unexpected text following translate"); }
1616 
1617 	else if(!strcmp(word,"copy")) {								// copy
1618 		CHECKS(fil,"need to enter filament name before copy");
1619 		itct=sscanf(line2,"%s",nm);
1620 		CHECKS(itct==1,"error reading filament name");
1621 		fil2=filAddFilament(fil->filtype,NULL,nm);
1622 		CHECKS(fil2,"failed to add filament");
1623 		er=filCopyFilament(fil2,fil);
1624 		CHECKS(!strnword(line2,2),"unexpected text following copy"); }
1625 
1626 	else {																				// unknown word
1627 		CHECKS(0,"syntax error within filament_type block: statement not recognized"); }
1628 
1629 	return fil;
1630 
1631 failure:
1632 	simParseError(sim,pfp);
1633 	return NULL; }
1634 
1635 
1636 /* filloadtype */
filloadtype(simptr sim,ParseFilePtr * pfpptr,char * line2)1637 int filloadtype(simptr sim,ParseFilePtr *pfpptr,char *line2) {
1638 	ParseFilePtr pfp;
1639 	char word[STRCHAR],errstring[STRCHAR];
1640 	int done,pfpcode,firstline2;
1641 	filamenttypeptr filtype;
1642 
1643 	pfp=*pfpptr;
1644 	done=0;
1645 	filtype=NULL;
1646 	firstline2=line2?1:0;
1647 	filenablefilaments(sim);
1648 
1649 	while(!done) {
1650 		if(pfp->lctr==0)
1651 			simLog(sim,2," Reading file: '%s'\n",pfp->fname);
1652 		if(firstline2) {
1653 			strcpy(word,"name");
1654 			pfpcode=1;
1655 			firstline2=0; }
1656 		else
1657 			pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
1658 
1659         *pfpptr=pfp;
1660         CHECKS(pfpcode!=3,"%s",errstring);
1661 
1662         if(pfpcode==0);																// already taken care of
1663         else if(pfpcode==2) {													// end reading
1664             done=1; }
1665         else if(!strcmp(word,"end_filament_type")) {		// end_filament_type
1666             CHECKS(!line2,"unexpected text following end_filament_type");
1667             return 0; }
1668         else if(!line2) {															// just word
1669             CHECKS(0,"unknown word or missing parameter"); }
1670         else {
1671             filtype=filtypereadstring(sim,pfp,filtype,word,line2);
1672             CHECK(filtype); }}													// failed but error has already been sent
1673 
1674 	CHECKS(0,"end of file encountered before end_filament_type statement");	// end of file
1675 
1676 	failure:																					// failure
1677 	if(ErrorType!=1) simParseError(sim,pfp);
1678 	*pfpptr=pfp=NULL;
1679 	return 1; }
1680 
1681 
1682 /* filloadfil */
filloadfil(simptr sim,ParseFilePtr * pfpptr,char * line2,filamenttypeptr filtype)1683 int filloadfil(simptr sim,ParseFilePtr *pfpptr,char *line2,filamenttypeptr filtype) {
1684 	ParseFilePtr pfp;
1685 	char word[STRCHAR],errstring[STRCHAR];
1686 	int done,pfpcode,firstline2;
1687 	filamentptr fil;
1688 
1689 	pfp=*pfpptr;
1690 	done=0;
1691 	firstline2=line2?1:0;
1692 	filenablefilaments(sim);
1693 	fil=NULL;
1694 
1695 	while(!done) {
1696 		if(pfp->lctr==0)
1697 			simLog(sim,2," Reading file: '%s'\n",pfp->fname);
1698 		if(firstline2) {
1699 			strcpy(word,"name");
1700 			pfpcode=1;
1701 			firstline2=0; }
1702 		else
1703 			pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
1704 
1705         *pfpptr=pfp;
1706         CHECKS(pfpcode!=3,"%s",errstring);
1707 
1708         if(pfpcode==0);																// already taken care of
1709         else if(pfpcode==2) {													// end reading
1710             done=1; }
1711         else if(!strcmp(word,"end_filament")) {				// end_filament_type
1712             CHECKS(!line2,"unexpected text following end_filament");
1713             return 0; }
1714         else if(!line2) {															// just word
1715             CHECKS(0,"unknown word or missing parameter"); }
1716         else {
1717             fil=filreadstring(sim,pfp,fil,filtype,word,line2);
1718             CHECK(fil); 															// failed but error has already been sent
1719             filtype=fil->filtype; }}
1720 
1721 	CHECKS(0,"end of file encountered before end_filament statement");	// end of file
1722 
1723 	failure:																					// failure
1724 	if(ErrorType!=1) simParseError(sim,pfp);
1725 	*pfpptr=pfp=NULL;
1726 	return 1; }
1727 
1728 
1729 /******************************************************************************/
1730 /******************************* structure updating ***************************/
1731 /******************************************************************************/
1732 
1733 
1734 /* filupdateparams */
filupdateparams(simptr sim)1735 int filupdateparams(simptr sim) {
1736 	(void) sim;
1737 	return 0; }
1738 
1739 
1740 /* filupdatelists */
filupdatelists(simptr sim)1741 int filupdatelists(simptr sim) {
1742 	(void) sim;
1743 	return 0; }
1744 
1745 
1746 /* filupdate */
filsupdate(simptr sim)1747 int filsupdate(simptr sim) {
1748 	int er;
1749 	filamentssptr filss;
1750 
1751 	filss=sim->filss;
1752 	if(filss) {
1753 		if(filss->condition<=SClists) {
1754 			er=filupdatelists(sim);
1755 			if(er) return er;
1756 			filsetcondition(filss,SCparams,1); }
1757 		if(filss->condition==SCparams) {
1758 			er=filupdateparams(sim);
1759 			if(er) return er;
1760 			filsetcondition(filss,SCok,1); }}
1761 	return 0; }
1762 
1763 
1764 /******************************************************************************/
1765 /**************************** core simulation functions ***********************/
1766 /******************************************************************************/
1767 
1768 
1769 /* filSegmentXSurface */
filSegmentXSurface(simptr sim,filamentptr fil,char endchar)1770 int filSegmentXSurface(simptr sim,filamentptr fil,char endchar) { //?? Currently only suitable for treadmilling, not dynamics
1771 	int s,p,mxs;
1772 	surfaceptr srf;
1773 	panelptr pnl;
1774 	double *pt1,*pt2,crosspt[3];
1775 	enum PanelShape ps;
1776 	filamenttypeptr filtype;
1777 	beadptr bead,beadplus;
1778 	segmentptr segment;
1779 
1780 	filtype = fil->filtype;
1781 
1782 	if(!sim->srfss) return 0;
1783 	if(endchar=='f') {
1784 		if(filtype->isbead){
1785 			bead=fil->beads[fil->frontbs];
1786 			beadplus=fil->beads[fil->frontbs+1];
1787 			pt1=bead->xyz;
1788 			pt2=beadplus->xyz; }
1789 		else{
1790 			segment=fil->segments[fil->frontbs];
1791 			pt1=segment->xyzfront;
1792 			pt2=segment->xyzback; }}
1793 	else {
1794 		if(filtype->isbead){
1795 			bead=fil->beads[fil->nbs+fil->frontbs-1];
1796 			beadplus=fil->beads[fil->nbs+fil->frontbs];
1797 			pt1=bead->xyz;
1798 			pt2=beadplus->xyz; }
1799 		else{
1800 			segment=fil->segments[fil->nbs+fil->frontbs];
1801 			pt1=segment->xyzfront;
1802 			pt2=segment->xyzback; }}
1803 
1804 	mxs=0;
1805 	for(s=0;s<sim->srfss->nsrf && !mxs;s++) {
1806 		srf=sim->srfss->srflist[s];
1807 		for(ps=(enum PanelShape)0;ps<PSMAX && !mxs;ps=(enum PanelShape)(ps+1))
1808 			for(p=0;p<srf->npanel[ps] && !mxs;p++) {
1809 				pnl=srf->panels[ps][p];
1810 				mxs=lineXpanel(pt1,pt2,pnl,3,crosspt,NULL,NULL,NULL,NULL,NULL,0); }}
1811 //				printf("pt1=(%g %g %g), pt2=(%g %g %g), mxs=%i\n",pt1[0],pt1[1],pt1[2],pt2[0],pt2[1],pt2[2],mxs); }}
1812 
1813 	return mxs; }
1814 
1815 
1816 /* filSegmentXFilament */
filSegmentXFilament(simptr sim,filamentptr fil,char endchar,filamentptr * filptr)1817 int filSegmentXFilament(simptr sim,filamentptr fil,char endchar,filamentptr *filptr) { //?? Currently only supported for segments not beads
1818 	int f,i,mxf,mn,mn1,ft;
1819 	double thk=0.0,*pt1,*pt2,dist=0.0;
1820 	filamentssptr filss;
1821 	filamenttypeptr filtype;
1822 	filamentptr fil2;
1823 	segmentptr segment;
1824 
1825 	if(endchar=='f') {
1826 		segment=fil->segments[fil->frontbs];
1827 		pt1=segment->xyzfront;
1828 		pt2=segment->xyzback;
1829 		thk=segment->thk;
1830 		mn=mn1=fil->frontbs;
1831 		if(fil->nbs>1) mn1=fil->frontbs+1; }
1832 	else {
1833 		segment=fil->segments[fil->nbs+fil->frontbs];
1834 		pt1=segment->xyzfront;
1835 		pt2=segment->xyzback;
1836 		mn=mn1=fil->nbs+fil->frontbs-1;
1837 		if(fil->nbs>1) mn1=fil->nbs+fil->frontbs-2; }
1838 
1839 	mxf=0;
1840 	filss=sim->filss;
1841 	fil2=NULL;
1842 	for(ft=0;ft<filss->ntype && !mxf;ft++) {
1843 		filtype=filss->filtypes[ft];
1844 		for(f=0;f<filtype->nfil && !mxf;f++) {
1845 			fil2=filtype->fillist[f];
1846 			for(i=fil2->frontbs;i<fil2->nbs+fil2->frontbs && !mxf;i++) {
1847 				if(!(fil2==fil && (i==mn || i==mn1))) {
1848 					dist=Geo_NearestSeg2SegDist(pt1,pt2,fil2->segments[i]->xyzfront,fil2->segments[i]->xyzfront);		//?? check
1849 					if(dist<thk+fil2->segments[i]->thk) mxf=1; }}}}
1850 	if(mxf && filptr)
1851 		*filptr=fil2;
1852 	return mxf; }
1853 
1854 
1855 /* filTreadmill */
filTreadmill(simptr sim,filamentptr fil,int steps)1856 void filTreadmill(simptr sim,filamentptr fil,int steps) {		//?? Currently for segments only
1857 	int i,mxs,done,iter;
1858 	double thk,length,angle[3],sigmamult,angletry[3],dcm[9];
1859 	filamentptr fil2;
1860 	filamenttypeptr filtype;
1861 
1862 	filtype=fil->filtype;
1863 	for(i=0;i<steps;i++) {
1864 		thk=fil->segments[0]->thk;
1865 		done=0;
1866 		sigmamult=1;
1867 		angletry[0]=angletry[1]=angletry[2]=0;
1868 		for(iter=0;iter<500 && !done;iter++) {
1869 			length=filRandomLength(filtype,thk,sigmamult);
1870 			filRandomAngle(filtype,thk,angle,sigmamult);
1871 			if(iter>0 && coinrandD(0.5))
1872 				filAddSegment(fil,NULL,length,angletry,thk,'b');
1873 			else
1874 				filAddSegment(fil,NULL,length,angle,thk,'b');
1875 			mxs=filSegmentXSurface(sim,fil,'b');
1876 			if(!mxs) {
1877 				mxs=filSegmentXFilament(sim,fil,'b',&fil2);
1878 				if(mxs) {
1879 					mxs=coinrandD(0.995);
1880 					Sph_DcmxDcmt(fil2->segments[fil2->nbs+fil2->frontbs-1]->adcm,fil->segments[fil->nbs+fil->frontbs-2]->adcm,dcm);
1881 					Sph_Dcm2Xyz(dcm,angletry); }}
1882 			if(mxs) {
1883 				filRemoveSegment(fil,'b');
1884 				sigmamult*=1.01; }
1885 			else done=1; }
1886 		if(done)
1887 			filRemoveSegment(fil,'f'); }
1888 
1889 	return; }
1890 
1891 
1892 /* filDynamics */
filDynamics(simptr sim)1893 int filDynamics(simptr sim) {
1894 	filamentssptr filss;
1895 	filamentptr fil;
1896 	filamenttypeptr filtype;
1897 	beadptr bead,beadplus,beadminus;
1898 	int f,b,d,ft;
1899 	double k1,k2; // Force & gaussian constants
1900 	int dim;
1901 
1902 	filss=sim->filss;
1903 	if(!filss) return 0;
1904 
1905 	dim=sim->dim;
1906 
1907 	for(ft=0;ft<filss->ntype;ft++) {
1908 		filtype=filss->filtypes[ft];
1909 		for(f=0;f<filtype->nfil;f++) {
1910 			fil=filtype->fillist[f];
1911 			if(filtype->treadrate>0)
1912 				filTreadmill(sim,fil,poisrandD(filtype->treadrate*sim->dt));
1913 
1914 			if(filtype->dynamics==FDrouse) {
1915 				k1 = 3*filtype->kT*sim->dt/(6*PI*filtype->viscosity*filtype->filradius*filtype->stdlen*filtype->stdlen); //?? Double check this is kuhn length squared
1916 				k2 = sqrt(2*filtype->kT/(6*PI*filtype->viscosity*filtype->filradius));
1917 
1918 				for(b=fil->frontbs;b<=fil->nbs+fil->frontbs;b++){
1919 					for(d=0;d<dim;d++){
1920 						fil->beads[b]->xyzold[d]=fil->beads[b]->xyz[d];}}
1921 					//?? PERHAPS include a check for number of segments to be >= 2
1922 
1923 				b=fil->frontbs;
1924 				bead=fil->beads[b];
1925 				beadplus=fil->beads[b+1];
1926 				for(d=0;d<dim;d++){
1927 					bead->xyz[d]=bead->xyzold[d]-k1*(bead->xyzold[d]-beadplus->xyzold[d])+k2*gaussrandD();}
1928 
1929 				for(b=fil->frontbs+1;b<fil->nbs+fil->frontbs;b++){
1930 					beadminus=fil->beads[b-1];
1931 					bead=fil->beads[b];
1932 					beadplus=fil->beads[b+1];
1933 					for(d=0;d<dim;d++){
1934 						bead->xyz[d]=bead->xyzold[d]-k1*(2*bead->xyzold[d]-beadminus->xyzold[d]-beadplus->xyzold[d])+k2*gaussrandD();}}
1935 
1936 				b=fil->nbs+fil->frontbs;
1937 				beadminus=fil->beads[b-1];
1938 				bead=fil->beads[b];
1939 				for(d=0;d<dim;d++){
1940 					bead->xyz[d]=bead->xyzold[d]-k1*(bead->xyzold[d]-beadminus->xyzold[d])+k2*gaussrandD();}
1941 				}}}
1942 
1943 	return 0; }
1944 
1945 
1946