Let's do a little (just a little) algebra.
Let x be the value in the central square; let x_i, i = 1, .., 8 index the values in the neighboring squares; and let r be the topographic ruggedness index. This recipe says r^2 equals the sum of (x_i - x)^2. Two things we can compute easily are (i) the sum of the values in the neighborhood, equal to s = Sum{ x_i } + x; and (ii) the sum of squares of the values, equal to t = Sum{ x_i^2 } + x^2. (These are focal statistics for the original grid and for its square.)
Expanding the squares gives
r^2 = Sum{ (x_i - x)^2 }
= Sum{ x_i^2 + x^2 - 2*x*x_i }
= Sum{ x_i^2 } + 8*x^2 - 2*x*Sum{x_i}
= [Sum{ x_i^2 } + x^2] + 7*x^2 - 2*x*[Sum{ x_i } + x - x]
= t + 7*x^2 - 2*x*[Sum{ x_i } + x] + 2*x^2
= t + 9*x^2 - 2*x*s.
For example, consider a neighborhood
1 2 3
4 5 6
7 8 9
Here, x = 5, s = 1+2+...+9 = 45, and t = 1+4+9+...+81 = 285. Then
(1-5)^2 + (2-5)^2 + ... + (9-5)^2 = 16 + 9 + 4 + 1 + 1 + 4 + 9 + 16 = 60 = r^2
and the algebraic equivalence says
60 = r^2 = 285 + 9*5^2 -2*5*45 = 285 + 225 - 450 = 60, which checks.
The workflow therefore is:
Given a DEM.
Compute s = Focal sum (over 3 x 3 square neighborhoods) of [DEM].
Compute DEM2 = [DEM]*[DEM].
Compute t = Focal sum (over 3 x 3 square neighborhoods) of [DEM2].
Compute r2 = [t] + 9*[DEM2] - 2*[DEM]*[s].
Return r = Sqrt([r2]).
This consists of 9 grid operations in toto, all of which are fast. They are readily carried out in the raster calculator (ArcGIS 9.3 and earlier), the command line (all versions), and Model Builder (all versions).
BTW, this is not an "average elevation change" (because elevation changes can be positive and negative): it is a root mean square elevation change. It is not equal to the "topographic position index" described at http://arcscripts.esri.com/details.asp?dbid=14156 , which (according to the documentation) equals x - (s - x)/8. In the example above, the TPI equals 5 - (45-5)/8 = 0 whereas the TRI, as we saw, is Sqrt(60).
I think I understand what you are trying to do now. You need to loop over all the meridional and zonal (24 each) TIFFs and use these as inputs for each raster in your directory.
If that is the case, then you need to nest iterators. The first one will loop over the raw images while the second loops over each M&Z file, passing this as input to the con function. You'll need an in-line variable to do that. The variable is the name of each M&Z file (which changes each loop).
If you want to loop over rasters in a directory, you need to use an iterator. You can run the model on all the rasters (restricted by image type if desired).
In simple, rough pseudo code, I would imagine it would look something like this:
rawpath = "directory containing raw images"
MZpath = "directory containing M&Z files"
for each raster in rawpath:
for each MZ in MZpath:
Con(raster,M,Z)
Con.save()
Edit: Apparently you can only have one iterator in modelbuilder. That's immensely annoying. As @Radar pointed out, you could nest models together. What I would do is get your model to work on just 1 image (using all the M&Z) files and then try to nest it within another model. Or learn Python, where this task would be very simple.
Best Answer
From what I know, it is not possible to use relative pixel index in ArcGIS raster calculator like you can do in GRASS. The best way to use such indexing method would be to use arcpy.rastertonumpyarray then loop on all pixel values to compute your angles. Something like
Of course, this code will not work as it is because it some indices will be out of range. Also, you should be careful because numpy indexing start at the top left.
Another workaround is to use a stack of shifted raster that you could use in raster calculator, but this would give you a lot of useless rasters.