trisurf-ng.git - Gitblit

Samo Penic

2016-05-15 fd8126e0965819f664ad0ae0c8d6a11e9509ea27

commit \| author \| age
7f6076	1	/* vim: set ts=4 sts=4 sw=4 noet : */
aec47d	2	#include<stdlib.h>
SP	3	#include<math.h>
	4	#include "general.h"
	5	#include "vertex.h"
	6	#include "bond.h"
	7	#include "triangle.h"
	8	#include "vesicle.h"
	9	#include "energy.h"
	10	#include "timestep.h"
	11	#include "cell.h"
	12	//#include "io.h"
9166cb	13	#include "io.h"
aec47d	14	#include<stdio.h>
SP	15	#include "vertexmove.h"
1ad6d1	16	#include <string.h>
43c042	17	#include "constvol.h"
aec47d	18
fedf2b	19	ts_bool single_verticle_timestep(ts_vesicle vesicle,ts_vertex vtx,ts_double *rn){
aec47d	20	ts_uint i;
SP	21	ts_double dist;
	22	ts_bool retval;
	23	ts_uint cellidx;
c0ae90	24	ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0;
ed31fe	25	ts_double costheta,sintheta,phi,r;
1ad6d1	26	//This will hold all the information of vtx and its neighbours
958e0e	27	ts_vertex backupvtx[20], constvol_vtx_moved=NULL, constvol_vtx_backup=NULL;
dcd350	28	memcpy((void )&backupvtx[0],(void )vtx,sizeof(ts_vertex));
a63f17	29
SP	30	//Some stupid tests for debugging cell occupation!
	31	/* cellidx=vertex_self_avoidance(vesicle, vtx);
	32	if(vesicle->clist->cell[cellidx]==vtx->cell){
	33	fprintf(stderr,"Idx match!\n");
	34	} else {
	35	fprintf(stderr,"***** Idx don't match!\n");
	36	fatal("ENding.",1);
	37	}
	38	*/
	39
352fad	40	//temporarly moving the vertex
672ae4	41	// vtx->x=vtx->x+vesicle->stepsize(2.0rn[0]-1.0);
SP	42	// vtx->y=vtx->y+vesicle->stepsize(2.0rn[1]-1.0);
	43	// vtx->z=vtx->z+vesicle->stepsize(2.0rn[2]-1.0);
	44
c0ae90	45	//random move in a sphere with radius stepsize:
ed31fe	46	r=vesicle->stepsize*rn[0];
M	47	phi=rn[1]2M_PI;
	48	costheta=2*rn[2]-1;
	49	sintheta=sqrt(1-pow(costheta,2));
672ae4	50	vtx->x=vtx->x+rsinthetacos(phi);
SP	51	vtx->y=vtx->y+rsinthetasin(phi);
	52	vtx->z=vtx->z+r*costheta;
	53
	54
c0ae90	55	//distance with neighbours check
8f6a69	56	for(i=0;i<vtx->neigh_no;i++){
352fad	57	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
8f6a69	58	if(dist<1.0 \|\| dist>vesicle->dmax) {
c0ae90	59	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
SP	60	return TS_FAIL;
8f6a69	61	}
aec47d	62	}
304510	63
M	64	// Distance with grafted poly-vertex check:
	65	if(vtx->grafted_poly!=NULL){
	66	dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
	67	if(dist<1.0 \|\| dist>vesicle->dmax) {
	68	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	69	return TS_FAIL;
	70	}
	71	}
	72
fe24d2	73	// TODO: Maybe faster if checks only nucleus-neighboring cells
M	74	// Nucleus penetration check:
	75	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	76	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	77	return TS_FAIL;
	78	}
	79
	80	//self avoidance check with distant vertices
	81	cellidx=vertex_self_avoidance(vesicle, vtx);
	82	//check occupation number
	83	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	84
aec47d	85	if(retval==TS_FAIL){
dcd350	86	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
aec47d	87	return TS_FAIL;
SP	88	}
1ad6d1	89
SP	90
c0ae90	91	//if all the tests are successful, then energy for vtx and neighbours is calculated
1ad6d1	92	for(i=0;i<vtx->neigh_no;i++){
dcd350	93	memcpy((void )&backupvtx[i+1],(void )vtx->neigh[i],sizeof(ts_vertex));
1ad6d1	94	}
aec47d	95
1121fa	96	if(vesicle->pswitch == 1 \|\| vesicle->tape->constvolswitch>0){
414b8a	97	for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
c0ae90	98	}
SP	99
	100	if(vesicle->tape->constareaswitch==2){
	101	for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area;
	102
	103	}
a63f17	104
aec47d	105	delta_energy=0;
fe5069	106
SP	107	// vesicle_volume(vesicle);
	108	// fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
43c042	109
aec47d	110	//update the normals of triangles that share bead i.
8f6a69	111	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17	112	oenergy=vtx->energy;
aec47d	113	energy_vertex(vtx);
a63f17	114	delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d	115	//the same is done for neighbouring vertices
8f6a69	116	for(i=0;i<vtx->neigh_no;i++){
SP	117	oenergy=vtx->neigh[i]->energy;
	118	energy_vertex(vtx->neigh[i]);
	119	delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
aec47d	120	}
414b8a	121
1121fa	122	if(vesicle->pswitch == 1 \|\| vesicle->tape->constvolswitch >0){
414b8a	123	for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
fbcbdf	124	if(vesicle->pswitch==1) delta_energy-=vesicle->pressure*dvol;
414b8a	125	};
43c042	126
c0ae90	127	if(vesicle->tape->constareaswitch==2){
SP	128	/* check whether the darea is gt epsarea */
	129	for(i=0;i<vtx->tristar_no;i++) darea+=vtx->tristar[i]->area;
	130	if(fabs(vesicle->area+darea-A0)>epsarea){
	131	//restore old state.
	132	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	133	for(i=0;i<vtx->neigh_no;i++){
	134	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
	135	}
	136	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
	137	//fprintf(stderr,"fajlam!\n");
	138	return TS_FAIL;
	139	}
	140
	141
	142	}
1121fa	143
SP	144	if(vesicle->tape->constvolswitch==2){
	145	/check whether the dvol is gt than epsvol /
	146	//fprintf(stderr,"DVOL=%1.16e\n",dvol);
	147	if(fabs(vesicle->volume+dvol-V0)>epsvol){
	148	//restore old state.
	149	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	150	for(i=0;i<vtx->neigh_no;i++){
	151	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
	152	}
	153	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
	154	//fprintf(stderr,"fajlam!\n");
	155	return TS_FAIL;
	156	}
	157
	158	} else
fe5069	159	// vesicle_volume(vesicle);
SP	160	// fprintf(stderr,"Volume before=%1.16e\n", vesicle->volume);
fbcbdf	161	if(vesicle->tape->constvolswitch == 1){
fe5069	162	retval=constvolume(vesicle, vtx, -dvol, &delta_energy_cv, &constvol_vtx_moved,&constvol_vtx_backup);
fbcbdf	163	if(retval==TS_FAIL){ // if we couldn't move the vertex to assure constant volume
SP	164	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	165	for(i=0;i<vtx->neigh_no;i++){
	166	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
	167	}
	168	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
fe5069	169	// fprintf(stderr,"fajlam!\n");
fbcbdf	170	return TS_FAIL;
SP	171	}
fe5069	172	// vesicle_volume(vesicle);
SP	173	// fprintf(stderr,"Volume after=%1.16e\n", vesicle->volume);
	174	// fprintf(stderr,"Volume after-dvol=%1.16e\n", vesicle->volume-dvol);
	175	// fprintf(stderr,"Denergy before=%e\n",delta_energy);
	176
fbcbdf	177	delta_energy+=delta_energy_cv;
fe5069	178	// fprintf(stderr,"Denergy after=%e\n",delta_energy);
fbcbdf	179	}
304510	180	/* No poly-bond energy for now!
fedf2b	181	if(vtx->grafted_poly!=NULL){
M	182	delta_energy+=
	183	(pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
	184	pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
	185	}
304510	186	*/
314f2d	187	// fprintf(stderr, "DE=%f\n",delta_energy);
aec47d	188	//MONTE CARLOOOOOOOO
SP	189	if(delta_energy>=0){
	190	#ifdef TS_DOUBLE_DOUBLE
3de289	191	if(exp(-delta_energy)< drand48())
aec47d	192	#endif
SP	193	#ifdef TS_DOUBLE_FLOAT
	194	if(expf(-delta_energy)< (ts_float)drand48())
	195	#endif
	196	#ifdef TS_DOUBLE_LONGDOUBLE
	197	if(expl(-delta_energy)< (ts_ldouble)drand48())
	198	#endif
	199	{
	200	//not accepted, reverting changes
fbcbdf	201	// fprintf(stderr,"MC failed\n");
dcd350	202	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
1ad6d1	203	for(i=0;i<vtx->neigh_no;i++){
a63f17	204	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1	205	}
SP	206
aec47d	207	//update the normals of triangles that share bead i.
dcd350	208	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
1ad6d1	209
fe5069	210	// fprintf(stderr, "before vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
43c042	211	if(vesicle->tape->constvolswitch == 1){
958e0e	212	constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
43c042	213	}
fe5069	214	// fprintf(stderr, "after vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
dd5aca	215	// vesicle_volume(vesicle);
SP	216	// fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
aec47d	217	return TS_FAIL;
SP	218	}
	219	}
2b14da	220	//accepted
fbcbdf	221	// fprintf(stderr,"MC accepted\n");
a63f17	222	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
SP	223	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	224	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	225	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	226	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
	227
	228	}
2b14da	229
1121fa	230	if(vesicle->tape->constvolswitch == 2){
SP	231	vesicle->volume+=dvol;
	232	} else
43c042	233	if(vesicle->tape->constvolswitch == 1){
fbcbdf	234	constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
43c042	235	}
c0ae90	236
SP	237	if(vesicle->tape->constareaswitch==2){
	238	vesicle->area+=darea;
	239	}
a63f17	240	// if(oldcellidx);
aec47d	241	//END MONTE CARLOOOOOOO
dd5aca	242	// vesicle_volume(vesicle);
SP	243	// fprintf(stderr,"Volume after success=%1.16e\n", vesicle->volume);
aec47d	244	return TS_SUCCESS;
SP	245	}
	246
fedf2b	247
M	248	ts_bool single_poly_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	249	ts_uint i;
	250	ts_bool retval;
	251	ts_uint cellidx;
304510	252	// ts_double delta_energy;
fedf2b	253	ts_double costheta,sintheta,phi,r;
304510	254	ts_double dist;
fedf2b	255	//This will hold all the information of vtx and its neighbours
M	256	ts_vertex backupvtx;
304510	257	// ts_bond backupbond[2];
fedf2b	258	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	259
	260	//random move in a sphere with radius stepsize:
	261	r=vesicle->stepsize*rn[0];
	262	phi=rn[1]2M_PI;
	263	costheta=2*rn[2]-1;
	264	sintheta=sqrt(1-pow(costheta,2));
	265	vtx->x=vtx->x+rsinthetacos(phi);
	266	vtx->y=vtx->y+rsinthetasin(phi);
	267	vtx->z=vtx->z+r*costheta;
	268
	269
	270	//distance with neighbours check
304510	271	for(i=0;i<vtx->neigh_no;i++){
M	272	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
	273	if(dist<1.0 \|\| dist>vesicle->dmax) {
	274	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	275	return TS_FAIL;
	276	}
	277	}
	278
	279	// Distance with grafted vesicle-vertex check:
	280	if(vtx==poly->vlist->vtx[0]){
	281	dist=vtx_distance_sq(vtx,poly->grafted_vtx);
	282	if(dist<1.0 \|\| dist>vesicle->dmax) {
	283	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	284	return TS_FAIL;
	285	}
	286	}
	287
fedf2b	288
M	289	//self avoidance check with distant vertices
	290	cellidx=vertex_self_avoidance(vesicle, vtx);
	291	//check occupation number
	292	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	293
	294	if(retval==TS_FAIL){
	295	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	296	return TS_FAIL;
	297	}
	298
	299
	300	//if all the tests are successful, then energy for vtx and neighbours is calculated
304510	301	/* Energy ignored for now!
fedf2b	302	delta_energy=0;
M	303	for(i=0;i<vtx->bond_no;i++){
	304	memcpy((void )&backupbond[i],(void )vtx->bond[i],sizeof(ts_bond));
	305
	306	vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
	307	bond_energy(vtx->bond[i],poly);
	308	delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
	309	}
	310
	311	if(vtx==poly->vlist->vtx[0]){
	312	delta_energy+=
	313	(pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
	314	pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
	315
	316	}
	317
	318
	319	if(delta_energy>=0){
	320	#ifdef TS_DOUBLE_DOUBLE
	321	if(exp(-delta_energy)< drand48() )
	322	#endif
	323	#ifdef TS_DOUBLE_FLOAT
	324	if(expf(-delta_energy)< (ts_float)drand48())
	325	#endif
	326	#ifdef TS_DOUBLE_LONGDOUBLE
	327	if(expl(-delta_energy)< (ts_ldouble)drand48())
	328	#endif
	329	{
	330	//not accepted, reverting changes
	331	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	332	for(i=0;i<vtx->bond_no;i++){
	333	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
	334	}
	335
	336	return TS_FAIL;
	337	}
	338	}
304510	339	*/
fedf2b	340
M	341	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
	342	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	343	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	344	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	345	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	346	}
	347	// if(oldcellidx);
	348	//END MONTE CARLOOOOOOO
	349	return TS_SUCCESS;
	350	}
58230a	351
M	352
	353
	354
	355	ts_bool single_filament_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	356	ts_uint i;
	357	ts_bool retval;
	358	ts_uint cellidx;
b30f45	359	ts_double delta_energy;
58230a	360	ts_double costheta,sintheta,phi,r;
M	361	ts_double dist[2];
	362	//This will hold all the information of vtx and its neighbours
b30f45	363	ts_vertex backupvtx,backupneigh[2];
58230a	364	ts_bond backupbond[2];
b30f45	365
M	366	//backup vertex:
58230a	367	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	368
	369	//random move in a sphere with radius stepsize:
	370	r=vesicle->stepsize*rn[0];
	371	phi=rn[1]2M_PI;
	372	costheta=2*rn[2]-1;
	373	sintheta=sqrt(1-pow(costheta,2));
	374	vtx->x=vtx->x+rsinthetacos(phi);
	375	vtx->y=vtx->y+rsinthetasin(phi);
	376	vtx->z=vtx->z+r*costheta;
	377
	378
	379	//distance with neighbours check
	380	for(i=0;i<vtx->bond_no;i++){
	381	dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
	382	if(dist[i]<1.0 \|\| dist[i]>vesicle->dmax) {
	383	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	384	return TS_FAIL;
	385	}
	386	}
	387
fe24d2	388	// TODO: Maybe faster if checks only nucleus-neighboring cells
M	389	// Nucleus penetration check:
	390	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	391	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	392	return TS_FAIL;
	393	}
	394
58230a	395
M	396	//self avoidance check with distant vertices
	397	cellidx=vertex_self_avoidance(vesicle, vtx);
	398	//check occupation number
	399	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	400	if(retval==TS_FAIL){
	401	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	402	return TS_FAIL;
	403	}
	404
	405	//backup bonds
	406	for(i=0;i<vtx->bond_no;i++){
	407	memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
	408	vtx->bond[i]->bond_length=sqrt(dist[i]);
	409	bond_vector(vtx->bond[i]);
b30f45	410	}
M	411
	412	//backup neighboring vertices:
	413	for(i=0;i<vtx->neigh_no;i++){
	414	memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a	415	}
M	416
	417	//if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45	418	delta_energy=0;
M	419
	420	if(vtx->bond_no == 2){
	421	vtx->energy = -(vtx->bond[0]->xvtx->bond[1]->x + vtx->bond[0]->yvtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length;
	422	delta_energy += vtx->energy - backupvtx.energy;
58230a	423	}
M	424
b30f45	425	for(i=0;i<vtx->neigh_no;i++){
M	426	if(vtx->neigh[i]->bond_no == 2){
	427	vtx->neigh[i]->energy = -(vtx->neigh[i]->bond[0]->xvtx->neigh[i]->bond[1]->x + vtx->neigh[i]->bond[0]->yvtx->neigh[i]->bond[1]->y + vtx->neigh[i]->bond[0]->z*vtx->neigh[i]->bond[1]->z)/vtx->neigh[i]->bond[0]->bond_length/vtx->neigh[i]->bond[1]->bond_length;
	428	delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
	429	}
58230a	430	}
M	431
b30f45	432	// poly->k is filament persistence length (in units l_min)
M	433	delta_energy *= poly->k;
58230a	434
M	435	if(delta_energy>=0){
	436	#ifdef TS_DOUBLE_DOUBLE
	437	if(exp(-delta_energy)< drand48() )
	438	#endif
	439	#ifdef TS_DOUBLE_FLOAT
	440	if(expf(-delta_energy)< (ts_float)drand48())
	441	#endif
	442	#ifdef TS_DOUBLE_LONGDOUBLE
	443	if(expl(-delta_energy)< (ts_ldouble)drand48())
	444	#endif
	445	{
	446	//not accepted, reverting changes
	447	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
b30f45	448	for(i=0;i<vtx->neigh_no;i++){
M	449	memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
	450	}
58230a	451	for(i=0;i<vtx->bond_no;i++){
b30f45	452	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
58230a	453	}
M	454
	455	return TS_FAIL;
	456	}
	457	}
	458
b30f45	459
58230a	460	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M	461	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	462	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	463	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	464	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	465	}
	466	// if(oldcellidx);
	467	//END MONTE CARLOOOOOOO
	468	return TS_SUCCESS;
	469	}