diff --git a/triangle_ccm_timeevol_plots.m b/triangle_ccm_timeevol_plots.m
new file mode 100644
index 0000000000000000000000000000000000000000..78d16e7455208a9a8bece590c91b4f29294a3f47
--- /dev/null
+++ b/triangle_ccm_timeevol_plots.m
@@ -0,0 +1,113 @@
+(* ::Package:: *)
+
+Needs["ErrorBarPlots`"]
+
+
+(* ::Section:: *)
+(*Plot successrate over time*)
+
+
+gsraw=Import[NotebookDirectory[]<>"data/graphdata_ccm_timeevol_canonical.m"];
+(* gsraw=SortBy[gsraw,{#[[1,1]]&,#[[1,2]]&}]; (* Sort by n and then by tau. The {} forces a *stable* sort because otherwise Mathematica sorts also on triangle count and other things. *) *)
+
+
+gdata=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}];
+(* Data format: *)
+(* gdata[[ tau index, n index, run index , datatype index ]] *)
+(* datatype index:
+1: {n,tau}
+2: edges
+3: HH tri seq
+4: {ccm1 #failed-attempts, ccm1 seq}
+5: {ccm2 #failed-attempts, ccm2 seq}
+*)
+nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];
+taulabels=Map["\[Tau] = "<>ToString[#]&,gdata[[All,1,1,1,2]]];
+
+
+(* ::Subsection:: *)
+(*Canonical dataset (old code for multiple ccm runs per run)*)
+
+
+(* For export *)
+selectedRun=gdata[[1,1,1]];
+measureSkip=1;
+numEdges=selectedRun[[2]];
+maxTris=Max[selectedRun[[3]]];
+maxTime=Length[selectedRun[[3]]];
+avgTris=Mean[selectedRun[[3,-maxTime/2;;-1]]];
+skipPts=Max[1,Round[maxTime/500]]; (* Plotting every point is slow. Plot only once per `skipPts` timesteps *)
+
+getCourseData[array_]:=Map[Mean,Partition[array,skipPts]];
+coarseData=getCourseData[selectedRun[[3]]];
+fullData=selectedRun[[3]];
+
+shiftedCCMdata[offset_,data_]:=If[Length[data]==0,Nothing,
+MapIndexed[{offset+(First[#2]-1)*skipPts,Mean[#1]}&,Partition[data,skipPts]]
+];
+coarseDatasCCM1=MapIndexed[shiftedCCMdata[numEdges*(First[#2]-1),#1]&,selectedRun[[4]]];
+coarseDatasCCM2=MapIndexed[shiftedCCMdata[numEdges*(First[#2]-1),#1]&,selectedRun[[5]]];
+
+getStartPoint[offset_,data_]:=If[Length[data]==0,
+{Red,Point[{offset,10}]},
+{Black,Point[{offset,data[[1]]}]}
+];
+startPoints1=MapIndexed[getStartPoint[numEdges*(First[#2]-1),#1]&,selectedRun[[4]]];
+
+plotlabel="n = "<>ToString[selectedRun[[1,1]]]<>", \[Tau] = "<>ToString[selectedRun[[1,2]]];
+
+plotTimeEvol1=ListPlot[coarseData,
+Joined->True,PlotRange->{0,3*avgTris},DataRange->{0,measureSkip*maxTime},
+Frame->True,FrameLabel->{"timesteps","number of triangles"},
+PlotLabel->plotlabel,ImageSize->300,
+PlotStyle->{Black},
+Epilog->{PointSize[0.015],startPoints1}
+];
+
+(*
+plotTimeEvol2=ListPlot[fullData,
+Joined->True,PlotRange\[Rule]All,DataRange->{0,measureSkip*maxTime},
+PlotStyle\[Rule]{Black,Opacity[0.3]}
+];*)
+
+plotTimeEvol3=ListPlot[coarseDatasCCM1,
+Joined->True,PlotRange->All
+];
+
+plotTimeEvol4=Show[plotTimeEvol1,plotTimeEvol3]
+(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)
+
+
+(* ::Subsection:: *)
+(*New code*)
+
+
+getCombinedData[run_]:=Module[{numEdges,maxTime,skipPts,getShiftedData,hhData,ccm1Data,ccm2Data},
+numEdges=run[[2]];
+maxTime=Length[run[[3]]];
+skipPts=Max[1,Round[maxTime/500]];
+
+getShiftedData[offset_,data_]:=If[Length[data]==0,{},
+MapIndexed[{offset+(First[#2]-1)*skipPts,Mean[#1]}&,Partition[data,skipPts]]
+];
+
+hhData=getShiftedData[0,run[[3]]];
+ccm1Data=getShiftedData[0*numEdges*run[[4,1]],run[[4,2]]];
+ccm2Data=getShiftedData[0*numEdges*run[[5,1]],run[[5,2]]];
+
+{Legended[hhData,"\[Tau] = "<>ToString[run[[1,2]]]],ccm1Data,ccm2Data}
+]
+
+dataSets=Map[getCombinedData,gdata,{3}];
+
+
+dataSetsFlattened=Flatten[dataSets,3];
+colorList=Table[ColorData[97,"ColorList"][[1+Floor[i/3]]],{i,0,Length[dataSetsFlattened]-1}];
+
+
+plot1=ListPlot[dataSetsFlattened,Joined->True,PlotRange->{{0,40000},{0,6000}},PlotStyle->colorList,ImageSize->300,PlotLabel->nlabels[[1]]]
+plot1log=ListLogPlot[dataSetsFlattened,Joined->True,PlotRange->{{0,40000},{1,15000}},PlotStyle->colorList,ImageSize->300,PlotLabel->nlabels[[1]]]
+
+
+Export[NotebookDirectory[]<>"plots/timeevol_ccm.pdf",plot1]
+Export[NotebookDirectory[]<>"plots/timeevol_ccm_log.pdf",plot1log]