After emailing the author, he described the process to generate this graph by saying:
"...we did simulations with these different SD values superimposed on the slopes and then computed power"
I've tried to replicate it but I'm having trouble getting similar results. So far what I've done is generated many linear models at each slope steepness (4xSD, 2xSD etc.) and for each number of sample size (3, 4 etc.) and then calculated the power of these regressions using the pwr.f2.test function in the pwr package in R.
For example, if I wanted to calculate the first square point on the graph (around (3,0.18)) I generated many sets of 3 points normally distributed around a line of a particular slope (determined by the SD chosen, 4xSD etc.). I then calculated the R^2 value of the fit of the points on the line and calculated the mean R^2 of all of the generated fits.
I then used the mean R^2 to calculate mean effect size (f^2 = R^2/(1-R^2)) and added it to the power equation:
pwr.2f.test(u, v, f2, sig) where u is the numerator degree of freedom (p - 1 = 2 - 1), v is the denominator df (n - p = 3 - 2 = 1), f2 is mean effect size and sig is level of significance (0.05) in this example.
Is this correct? The power values I get, however tend to 1 much quicker than this graph. I've also tried calculating the power for each model generated individually and getting the mean of all the powers but these results seem even further away from the graph. Sorry if this is incoherent, I'm relatively new to statistics.