diff --git a/triangle_ecm_initialtris.m b/triangle_ecm_initialtris.m
new file mode 100644
index 0000000000000000000000000000000000000000..f7f5382593f7840d05745fd0f45cc99f3109839f
--- /dev/null
+++ b/triangle_ecm_initialtris.m
@@ -0,0 +1,204 @@
+(* ::Package:: *)
+
+Quit[]
+
+
+Needs["ErrorBarPlots`"]
+
+
+(* ::Section:: *)
+(*Data import*)
+
+
+gsraw=Import[NotebookDirectory[]<>"data/graphdata_ecm_initialtris.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]]&}];
+(* Data format: *)
+(* gdata[[ tau index, n index, datatype index ]] *)
+(* datatype index:
+1: {n,tau}
+2: {uniform triangle samples}
+3: {ECM triangle samples}
+*)
+
+
+(* ::Section:: *)
+(*Erased configuration model*)
+
+
+(* ::Subsection:: *)
+(*Distribution of initial #triangles for ECM compared to uniform triangle distribution*)
+
+
+getHistogram[run_]:=Histogram[{run[[2]],run[[3]]},Automatic,"Probability",
+ChartLegends->Placed[{"Uniform","ECM"},Bottom],
+ImageSize->250,
+Frame->True,
+FrameLabel->{"Triangles","Probability"},
+PlotLabel->("n = "<>ToString[run[[1,1]]]<>", \[Tau] = "<>ToString[run[[1,2]]])
+];
+histograms=Map[getHistogram,gdata,{2}];
+TableForm[histograms]
+
+
+(* ::Subsubsection:: *)
+(*Exporting plots*)
+
+
+Needs["MaTeX`"]
+getHistogram[run_,bins_,plotrange_,paddings_,tickDelta_,bottomTicks_,textpos_,framelabel_]:=Histogram[{run[[3]],run[[2]]},bins,"Probability",
+ImageSize->150+paddings[[1,1]]+paddings[[1,2]],
+ImagePadding->paddings,
+AspectRatio->4/6,
+PlotRange->plotrange,
+Frame->True,
+FrameLabel->framelabel,
+FrameTicks->{{{#,NumberForm[#,{2,2}]}&/@Range[0,1,tickDelta],Automatic},{bottomTicks,Automatic}},
+Epilog->Text[MaTeX["n = "<>ToString[run[[1,1]]]<>",\\; \\tau = "<>ToString[run[[1,2]]]],textpos],
+ChartStyle->{ColorData[97][1],ColorData[97][2]}
+];
+textpos=Scaled[{0.65,0.90}];
+ticks1={#,ToString[#]}&/@Range[800,5000,800];
+ticks2={#*10^4,ToString[#\[CenterDot]Superscript[10,4],TraditionalForm]}&/@Range[2,10,2];
+{h1,h2,h3}={
+getHistogram[gdata[[2,1]],{50},{0,0.25},{{40,0},{12,0}},0.10,ticks1,textpos,{None,"Probability"}],
+getHistogram[gdata[[6,1]], {3},{0,0.10},{{40,0},{12,0}},0.05,Automatic,textpos,{None,"Probability"}],
+getHistogram[gdata[[10,1]],{1},{0,0.18},{{40,0},{32,0}},0.05,Automatic,textpos,{"Triangles","Probability"}]
+};
+(*plotgrid1=Column[{h1,h2,h3}]*)
+{h4,h5,h6}={
+getHistogram[gdata[[2,-1]],{800},{0,0.50},{{20,5},{12,0}},0.10,ticks2,textpos,None],
+getHistogram[gdata[[6,-1]], {10},{0,0.10},{{20,5},{12,0}},0.05,Automatic,textpos,None],
+getHistogram[gdata[[10,-1]], {1},{0,0.10},{{20,5},{32,0}},0.05,Automatic,textpos,{"Triangles"}]
+};
+(*SwatchLegend[{ColorData[97][1],ColorData[97][2]},{"ECM","Uniform"},LegendLayout\[Rule]"Row"]*)
+plotgrid2=Grid[Transpose[{{SwatchLegend[{ColorData[97][1]},{"ECM"}],h1,h2,h3},{SwatchLegend[{ColorData[97][2]},{"Uniform"}],h4,h5,h6}}]]
+
+
+Export[NotebookDirectory[]<>"plots/ecm_initialtris2.pdf",plotgrid2]
+
+
+(* ::Section:: *)
+(*As function of n*)
+
+
+dataPointsUniform=Map[{#[[1,1]],Mean[#[[2]]]}&,gdata,{2}];
+dataPointsECM=Map[{#[[1,1]],Mean[#[[3]]]}&,gdata,{2}];
+taulabels=Map["\[Tau] = "<>ToString[#[[1,1,2]]]&,gdata];
+
+(* Standard Deviation with division by N instead of N-1 *)
+mySD[xs_]:=Sqrt[Total[(xs-Mean[xs])^2]/Length[xs]];
+(* { {x,y}, ErrorBar[err] } *)
+getErrorBars[run_]:=Module[{n,tau,avgUni,avgECM,sdUni,sdECM},
+{n,tau}=run[[1]];
+avgUni=Mean[run[[2]]];
+sdUni=mySD[run[[2]]];
+avgECM=Mean[run[[3]]];
+sdECM=mySD[run[[3]]];
+{
+{{n,avgUni},ErrorBar[sdUni]},
+{{n,avgECM},ErrorBar[sdECM]}
+}
+]
+allErrorBars=Map[getErrorBars,gdata,{2}];
+
+
+Map[ErrorListPlot[Transpose[#],
+ImageSize->500,
+Frame->True,
+PlotMarkers->Automatic
+]&,allErrorBars[[{1,6,11}]]]
+
+
+(* ::Subsection:: *)
+(*Fitting the log-log-plot*)
+
+
+nRange=All;
+
+(* Weight: 1/err^2 *)
+(* Weight: N/Total[(xs-Mean[xs])^2] *)
+getWeight[xs_]:=1/(Log[Mean[xs]]-Log[Mean[xs]-Sqrt[Total[(xs-Mean[xs])^2]/Length[xs]]])^2;
+(* Several runs for fixed tau but different n *)
+getFitData[runs_,index_]:=Map[{Log[#[[1,1]]],Log[Mean[#[[index]]]],getWeight[#[[index]]]}&,runs];
+
+uniformFitData=Map[getFitData[#,2]&,gdata[[All,nRange]]];
+ECMFitData=Map[getFitData[#,3]&,gdata[[All,nRange]]];
+
+uniformFits=Map[LinearModelFit[#[[All,{1,2}]],logn,logn(*,Weights\[Rule]#[[All,3]]*)]&,uniformFitData];
+ECMFits=    Map[LinearModelFit[#[[All,{1,2}]],logn,logn(*,Weights\[Rule]#[[All,3]]*)]&,ECMFitData];
+(*
+uniformloglog=Log[dataPointsUniform[[All,nRange]]];
+ECMloglog=Log[dataPointsECM[[All,nRange]]];
+
+uniformFits=Map[LinearModelFit[#,logn,logn]&,uniformloglog];
+ECMFits=Map[LinearModelFit[#,logn,logn]&,ECMloglog];
+*)
+uniformFuncs=Map[#[logn]&,uniformFits];
+ECMFuncs=Map[#[logn]&,ECMFits];
+
+
+uniformFits[[1]]["ParameterTable"] (* Get `Standard Error' by "ParameterErrors" *)
+uniformFits[[1]]["ParameterConfidenceIntervalTable"] (* Get confidence by "ParameterConfidenceIntervals *)
+(* estimate +- standard error *)
+uniformFits[[1]]["BestFitParameters"]-uniformFits[[1]]["ParameterErrors"]
+uniformFits[[1]]["BestFitParameters"]+uniformFits[[1]]["ParameterErrors"]
+
+
+tauChoices={1,4,6,8,11};
+taulabels=Map["\[Tau] = "<>ToString[NumberForm[#[[1,1,2]],{3,2}]]&,gdata[[tauChoices]]];
+
+repeatColors[n_,k_]:=Table[ColorData[97][Mod[i,k,1]],{i,1,n}]
+
+plot1=Show[ListLogLogPlot[Evaluate[dataPointsUniform[[tauChoices]]~Join~dataPointsECM[[tauChoices]]],
+Frame->True,
+FrameTicks->{{Table[{10^k,Superscript[10,k]},{k,0,6}],Automatic},{{1000,2000,5000,10000},Automatic}},
+FrameLabel->{"n","triangles"},
+ImageSize->250,
+PlotMarkers->Automatic,
+PlotLegends->taulabels,
+PlotStyle->repeatColors[10,5]
+],
+Plot[Evaluate[uniformFuncs[[tauChoices]]],{logn,1,20000}],
+Plot[Evaluate[ECMFuncs[[tauChoices]]],{logn,1,20000}]
+]
+
+
+Export[NotebookDirectory[]<>"plots/avgtris_n.pdf",plot1]
+
+
+(* ::Subsection:: *)
+(*T(\[Tau]) including error bars*)
+
+
+gsraw2=Import[NotebookDirectory[]<>"data/graphdata_exponent_hightau.m"];
+gsraw2=SortBy[gsraw2,#[[1,1]]&]; (* Sort by n *)
+averagesGrouped=GatherBy[gsraw2,{#[[1,2]]&,#[[1,1]]&}];
+averagesLoglogdata=Map[{Log[#[[1,1,1]]],Log[Mean[#[[All,2]]]]}&,averagesGrouped[[All,nRange]],{2}];
+averagesFitsExtra=Map[LinearModelFit[#,logn,logn]&,averagesLoglogdata];
+avgTauValues=averagesGrouped[[All,1,1,1,2]];
+averagesExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,
+Transpose[{avgTauValues-0.000,averagesFitsExtra}]];
+
+
+tauValues=gdata[[All,1,1,2]];
+
+(* For visual, shift the tau values slightly left or right to distinguish the two datasets *)
+uniformExponents=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&, Transpose[{tauValues+0.000,uniformFits}]];
+ECMExponents    =Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&, Transpose[{tauValues+0.000,ECMFits}]];
+
+Needs["MaTeX`"]
+
+plot2=Show[
+ErrorListPlot[{ECMExponents,uniformExponents,averagesExponentsErrorBars},
+Joined->True,PlotMarkers->Automatic,
+PlotLegends->Placed[{"ECM","canonical","average"},{Left,Bottom}],
+Frame->True,FrameLabel->{"tau","triangle powerlaw exponent"},
+PlotRange->{{2,3},{0,1.6}},
+ImageSize->300],
+Plot[3/2(3-tau),{tau,2,3},PlotStyle->{Black,Dashed},PlotLegends->Placed[LineLegend[{MaTeX["\\frac{3}{2}(3-\\tau)"]},LegendMarkerSize->20],{Left,Bottom}]]]
+
+
+Export[NotebookDirectory[]<>"plots/triangle_exponent.pdf",plot2]