(* Clean version for generating data for paper *)

Definitions

Derivatives, Curvature and Forces

(* d1, d2, d3 and d4 are derivatives of a discretised vector function r *)
d1[r_]:=(Length[r]-1)Join@@{{-(1/2)r[[3]]+2r[[2]]-(3/2)r[[1]]},((Differences[r,1,2])/2),{(1/2)r[[-3]]-2r[[-2]]+(3/2)r[[-1]]}};

d2[r_]:=(((Length[r]-1))^2)Join@@{{-r[[4]]+4r[[3]]-5r[[2]]+2r[[1]]},Differences[r,2],{-r[[-4]]+4r[[-3]]-5r[[-2]]+2r[[-1]]}};

d3[r_]:=(((Length[r]-1))^3)Join@@{{-(3/2)r[[5]]+7r[[4]]-12r[[3]]+9r[[2]]-(5/2)r[[1]],-(3/2)r[[6]]+7r[[5]]-12r[[4]]+9r[[3]]-(5/2)r[[2]]},
((Differences[r,1,4])/2)-Rest@Most@Differences[r,1,2],
{(3/2)r[[-6]]-7r[[-5]]+12r[[-4]]-9r[[-3]]+(5/2)r[[-2]],(3/2)r[[-5]]-7r[[-4]]+12r[[-3]]-9r[[-2]]+(5/2)r[[-1]]}}; 

d4[r_]:=(((Length[r]-1))^4)Join@@{{-2r[[6]]+11r[[5]]-24r[[4]]+26r[[3]]-14r[[2]]+3r[[1]],-2r[[7]]+11r[[6]]-24r[[5]]+26r[[4]]-14r[[3]]+3r[[2]]},Differences[r,4],{-2r[[-7]]+11r[[-6]]-24r[[-5]]+26r[[-4]]-14r[[-3]]+3r[[-2]],-2r[[-6]]+11r[[-5]]-24r[[-4]]+26r[[-3]]-14r[[-2]]+3r[[-1]]}};

(* Magnitude of d1 *)
d1Mag[r_]:=(Length[r]-1)Sqrt@MapThread[Dot,{((d1[r])/(Length[r]-1)),((d1[r])/(Length[r]-1))}];

(* Unit tangent and normal vectors *)
That[r_]:=d1[r]/d1Mag[r];
Nhat[r_]:=RotationTransform[π/2]/@That[r];

(* Derivative of the magnitude of d1 *)
d1Magd[r_]:=MapThread[Dot,{d2[r],That[r]}];

(* Curvature and its derivatives *)
k[r_]:=MapThread[Dot,{((d2[r])/(d1Mag[r])),Nhat[r]}];
kd[r_]:=MapThread[Dot,{((d3[r])/(d1Mag[r])),Nhat[r]}]-2MapThread[Dot,{((d2[r])/(d1Mag[r])),That[r]}]MapThread[Dot,{((d2[r])/(d1Mag[r])),Nhat[r]}];
kdd[r_]:=MapThread[Dot,{((d4[r])/(d1Mag[r])),Nhat[r]}]-3(MapThread[Dot,{((d3[r])/(d1Mag[r])),That[r]}]MapThread[Dot,{((d2[r])/(d1Mag[r])),Nhat[r]}]+MapThread[Dot,{((d2[r])/(d1Mag[r])),That[r]}]MapThread[Dot,{((d3[r])/(d1Mag[r])),Nhat[r]}])-2MapThread[Dot,{((d2[r])/(d1Mag[r])),Nhat[r]}](((MapThread[Dot,{((d2[r])/(d1Mag[r])),Nhat[r]}])^2)-3((MapThread[Dot,{((d2[r])/(d1Mag[r])),That[r]}])^2));

(* Stretching forces *)
FSbulk[r_,gS_]:=gS(d2[r]-k[r]Nhat[r]);
FSend[r_,gS_]:=(Length[r]-1)gS(d1Mag[r]-1)That[r];
FS[r_,gS_]:=ReplacePart[FSbulk[r,gS],
{1->FSend[r,gS][[1]],-1->-FSend[r,gS][[-1]]}];

(* Bending forces *)
FBbulk[r_,gB_]:=gB(kd[r](MapThread[Dot,{((d2[r])/((d1Mag[r])^2)),Nhat[r]}]That[r]+MapThread[Dot,{((d2[r])/((d1Mag[r])^2)),That[r]}]Nhat[r])-((kdd[r])/(d1Mag[r]))Nhat[r]);
FBend[r_,gB_]:=(((Length[r]-1)gB)/(d1Mag[r]))(((k[r])^2)That[r]+(k[r]MapThread[Dot,{((d2[r])/(d1Mag[r])),That[r]}]-2kd[r])Nhat[r]);
FB[r_,gB_]:=ReplacePart[FBbulk[r,gB],
{1->FBend[r,gB][[1]],-1->-FBend[r,gB][[-1]]}];

(* Active force *)
FA[r_,α_]:=Cos[α]That[r]+Sin[α]Nhat[r];

(* Total force *)
F[r_,gS_,gB_,α_]:=FS[r,gS]+FB[r,gB]+FA[r,α];

Generic Simulation Functions

(* Step forward using a Backward Euler Method, time step δt,
wp = working precision to use for calculations *)
stepBE[r_,gS_,gB_,α_,δt_,addNoiseQ_:False,wp_:MachinePrecision]:=
Module[{equs,vars,n=Length[r]},

equs=Flatten[(#-r-F[#,gS,gB,α]δt)&@Table[{x[i],y[i]},{i,1,n}]];
If[addNoiseQ,
vars=Flatten[Table[{{x[i],r[[i,1]]+RandomReal[(({-1,1})/(n^3)),WorkingPrecision->wp]},{y[i],r[[i,2]]+RandomReal[(({-1,1})/(n^3)),WorkingPrecision->wp]}},{i,1,n}],1],
vars=Flatten[Table[{{x[i],r[[i,1]]},{y[i],r[[i,2]]}},{i,1,n}],1]
];

Table[{x[i],y[i]},{i,1,n}]/.FindRoot[equs,vars,WorkingPrecision->wp]
]


(* An algorithm which uses an addaptive time step to compute the evolution of the curve.
tMax = maximum time to evolve over
nStepsMax = maximum number of steps to try. 
Failed steps count!
Ends when either condition is met.
δtMax = maximum time step to ever try. If leftout, default to (10^(-2)).
δtMin = min. time step to try. If a smaller time step is "requested", simulation immediately ends. Default to (10^(-13)).
wp = working precision to use for computations.
If the filament tends to a constant shape (as measured by the spatial mean of (∂_t)k), the simulation can end after a few simulation seconds. δk and tConstMax determine the max change allowed to be considered constant and the amount of time that shoult pass with constant shape
 *)
addaptiveBE[r0_,gS_,gB_,α_,tMax_,nStepsMax_,δtMax_:(10^(-2)),δtMin_:(10^(-13)),wp_:MachinePrecision,δk_:(10^(-3)),tConstMax_:2]:=
First@Module[{δt=δtMax,t=0,stepN=0,
r1=r0,r2,addNoiseQ=False,
k1=k[r0],k2,tConst=0
},
Reap[
Sow[{t,r0,k1,d1Mag[r0]}];
While[t<tMax&&stepN<nStepsMax&&δt>=δtMin&&tConst<tConstMax,
stepN++;
r2=Quiet@Check[stepBE[r1,gS,gB,α,δt,addNoiseQ,wp],0,{FindRoot::cvmit,FindRoot::lstol}];

If[r2===0||Total[Norm/@(r2-r1)]<δt,
δt/=10,
t+=δt;k2=k[r2];
Sow[{t,r2,k2,d1Mag[r2]}];
If[Mean@Abs[((k2-k1)/(δt))]<=δk,tConst+=δt,tConst=0];
r1=r2;k1=k2;If[δt<δtMax,δt*=10];
];
]][[2]]
]

(*addaptiveBE[r0_,gS_,gB_,α_,tMax_,nStepsMax_,δtMax_:(10^(-2)),δtMin_:(10^(-13)),wp_:MachinePrecision]:=
First@Module[{δt=δtMax,t=0,stepN=0,
r1=r0,r2,addNoiseQ=False
},
Reap[
Sow[{t,r0}];
While[t<tMax&&stepN<nStepsMax&&δt>=δtMin,
stepN++;
r2=Quiet@Check[stepBE[r1,gS,gB,α,δt,addNoiseQ,wp],0,{FindRoot::cvmit,FindRoot::lstol}];

If[r2===0||Total[Norm/@(r2-r1)]<δt,
δt/=10,
t+=δt;Sow[{t,r2}];r1=r2;
If[δt<δtMax,δt*=10]
];


]][[2]]
]*)

Initial Condition Functions

(* Functions to return predicted w0 given gS, gB and α. If no solution exists, returns w0 = 0 *)
w0Large[gS_,gB_,α_]:=
With[{roots=SolveValues[(gB/(Sin[α]))(w^4)==w-((Sin[α])/gS),w,Reals]},
If[Length[roots]==0,0,Max[roots]]
];
w0Small[gS_,gB_,α_]:=
With[{roots=SolveValues[(gB/(Sin[α]))(w^4)==w-((Sin[α])/gS),w,Reals]},
If[Length[roots]==0,0,Min[roots]]
];

(* Functions to return predicted u0 given gS, gB and α. If no w0 solution exists, returns u0 = 1 *)
u0Large[gS_,gB_,α_]:=
With[{w0=w0Large[gS,gB,α]},
If[w0==0,1,(gB/(Sin[α]))(w0^3)]];
u0Small[gS_,gB_,α_]:=
With[{w0=w0Small[gS,gB,α]},
If[w0==0,1,(gB/(Sin[α]))(w0^3)]];

(* Create uniformly stretched and uniformly curved initial condition with nPoints.
Passes through origin *)
constCurve[w0_,u0_,nPoints_]:=
If[w0==0,
Subdivide[{0,-(u0/2)},{0,(u0/2)},nPoints-1],
(({Cos[#]-1,Sin[#]})/w0)&/@Subdivide[-((w0 u0)/2),((w0 u0)/2),nPoints-1]
];

Specific Simulation Functions

(* Auxillary Function for creating filenames.
Takes a number n and returns a string mEk where m and k are integers. m is returned with p digits (p precision) To extract n, use n = m*10^k  *)
numString[n_,p_]:=ToString[Round[First[MantissaExponent[n]]*(10^p)]]<>"E"<>ToString[Last[MantissaExponent[n]]-p];

(* Create data for constant curvature initial condition *)
createDataConstCurvatureLarge[gS_,gB_,α_,nPoints_,tMax_,nStepsMax_,δtMax_:(10^(-2)),δtMin_:(10^(-13)),wp_:MachinePrecision,δk_:(10^(-3)),tConstMax_:2]:=
Module[
{r0,data,assoc,filename},
r0=N[constCurve[w0Large[gS,gB,α],u0Large[gS,gB,α],nPoints],wp];
data=addaptiveBE[r0,gS,gB,α,tMax,nStepsMax,δtMax,δtMin,wp,δk,tConstMax];
assoc=<| 
"times"-> data[[All,1]],
"filaments"->data[[All,2]],
"kdata"->data[[All,3]],
"drdata"->data[[All,4]],
"gS"->gS,
"gB"->gB,
"alpha"->α
|>;
filename="constCurvLarge_gS_"<>numString[gS,3]<>"_gB_"<>numString[gB,3]<>"_a_"<>numString[α,1]<>".mx";
Export[filename,assoc,OverwriteTarget->"KeepBoth"];
];

Data Generation

SetDirectory[NotebookDirectory[]];

αVal=(10^(-1));
gSVal=(10^1)Sin[αVal];
gBVals=2^Range[-13+10/23,0,10/23]Sin[αVal];
gBVals=Table[(2^(-13-i(10/23))),{i,5,8}]Sin[αVal];

nPointsVal=250;
tMaxVal=8;
nStepsMaxVal=8000; 
(* Plot showing initial curve states *)
Table[
constCurve[w0Large[gSVal,gBVal,αVal],u0Large[gSVal,gBVal,αVal],nPointsVal],{gBVal,gBVals}];
ListPlot[%,PlotRange->{{-1,1},{-1,1}},AspectRatio->Automatic,Joined->True]

(* Generate data in parallel *)
δtMaxVal=(10^(-2));
δtMinVal=(10^(-8));
ParallelTable[createDataConstCurvatureLarge[gSVal,gBVal,αVal,nPointsVal,tMaxVal,nStepsMaxVal,δtMaxVal,δtMinVal],{gBVal,gBVals}];

(* α = 0.1 , gS/Sin[α] = 10, gB/Sin[α] variable *)
αVal=0.1;
gSVal=(10^1)Sin[αVal];
gBVals=10^Subdivide[Log10[2 (10^(-3))],Log10[((3^3)/(4^4))((((gSVal/(Sin[αVal]))))^3)],8-1]Sin[αVal];
nPointsVal=250;
tMaxVal=5;
nStepsMaxVal=2500; 
(* nPointsVal = 250 
& 
wp = MachinePrecision -> ~30 sec/step ?*)

(* Plot showing initial curve states *)
Table[
constCurve[w0Large[gSVal,gBVal,αVal],u0Large[gSVal,gBVal,αVal],nPointsVal],{gBVal,gBVals}];
ListPlot[%,PlotRange->{{-1,1},{-1,1}},AspectRatio->Automatic,Joined->True]

(* Generate data in parallel *)
ParallelTable[createDataConstCurvatureLarge[gSVal,gBVal,αVal,nPointsVal,tMaxVal,nStepsMaxVal],{gBVal,Reverse@gBVals}];

(* α = 0.1 , gS/Sin[α] = (10^2), gB/Sin[α] variable *)
αVal=0.1;
gSVal=(10^2)Sin[αVal];
gBVals=10^Subdivide[Log10[2 (10^(-3))],Log10[((3^3)/(4^4))((((gSVal/(Sin[αVal]))))^3)],8-1]Sin[αVal];
nPointsVal=250;
tMaxVal=5;
nStepsMaxVal=2500; 
(* nPointsVal = 250 
& 
wp = MachinePrecision -> ~30 sec/step ?*)

(* Plot showing initial curve states *)
Table[
constCurve[w0Large[gSVal,gBVal,αVal],u0Large[gSVal,gBVal,αVal],nPointsVal],{gBVal,gBVals}];
ListPlot[%,PlotRange->{{-1,1},{-1,1}},AspectRatio->Automatic,Joined->True]

(* Generate data in parallel *)
ParallelTable[createDataConstCurvatureLarge[gSVal,gBVal,αVal,nPointsVal,tMaxVal,nStepsMaxVal],{gBVal,Reverse@gBVals}];

Available Data & Additional Processing

SetDirectory[NotebookDirectory[]<>"Old data/13"];
fileList=FileNames["*.mx"];

Visualise Data

(* visualise some data *)
data=N@Import[fileList[[11]]];

(* simulation parameters *)
TableForm[{data["alpha"],data["gS"]/Sin[data["alpha"]],data["gB"]/Sin[data["alpha"]],data["times"][[-1]],Length[data["times"]]},TableHeadings->{{"α","gS/Sin(α)","gB/Sin(α)","tMax","nSteps"},None}]

(* Filament evolution *)
Manipulate[
ListPlot[{data["filaments"][[1]],data["filaments"][[i]]},Joined->True,PlotRange->2{{-1,1},{-1,1}},AspectRatio->Automatic,Epilog->{Text[ToString[N@data["times"][[i]]],Scaled@{0.05,0.95}]}],
{i,1,Length[data["filaments"]],1,Appearance->"Open"}]

(* Curvature evolution *)
Manipulate[
ListPlot[data["kdata"][[i]],
Joined->True,PlotRange->{0,4},
Epilog->{Text[ToString[N@data["times"][[i]]],Scaled@{0.05,0.95}]}],
{i,1,Length[data["kdata"]],Appearance->"Open"}]