Loop over orbital elements

gistfile1.r

# Create grids for periods, masses, ecc, mean anomaly, and longitude of pericenter
periods <- ...
masses <- ...
eccentricities <- c(0, 0.1, 0.2)
mas <- seq(0, 360, length.out=(N/2))[-1]
lops <- seq(0, 360, length.out=(N/2))[-1]

# Create a list of detection matrices; each element corresponds to the recovery rate matrix
# in (P, M) for the eccentricity.
detection.list <- list()

for e in eccentricities {
    # If eccentricity is circular, then only one grid point is needed for the l.o.p. (0),
    # otherwise use the one defined at the top of the file
    if (e == 0)
        loop.lops <- 0
    else
        loop.lops <- lops

    # Create the (P, M) matrix
    f <- matrix(nrow = N, ncol = N)
    
    # loop over periods and masses, as
    for P in periods {
        # ... skip if P > detectable.P, as usual ...
        
        for M in masses {
            
            detections <- 0
            # loop over mas and loop.lops (which is either 0 or the lop vector)
            for MA in mas {                
                for w in loop.lops {
                    # [ create dataset corresponding to signal due to planet with (P, M, MA, e, w) ] 
                    if (detectable)
                        detections <- detections + 1                    
                }
            }
            ...

            f[i, j] <- detections / (length(mas) * length(loop.lops))

        }
    }

    # append the matrix to the list
    detection.list[[length(detection.list)+1]] <- f
}