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