Since I've been fooling around with R lately, that's what I used. My code is below. It is a very simplified model, simulating one flat disk transiting another. The code produces random points in a box around the "planet" and discards the ones not within the planet's radius. Those points are then moved to the appropriate place in the transit across the "star", as determined by the step in the loop. The points are checked to see if they lay within the star's radius, and that count is divided by the total number of points to determine the proportion of the planet's area that is intersecting the star. That area is then divided by the star's area to find how much of the star's light is blocked. Those proportions are collected at each point of the planet's trip and can then a simulated light curve can be charted.
The ratio of planet to star radius can be changed, as can the crossing latitude. The simulation parameters of number of points generated and transit step length can also be changed. Again, this is nothing fancy and there are probably much better ways of doing this, but it was a fun project to chew on.
star.radius <- 1 # Conceptually, it's easier if this is 1
planet.radius <- .20
star.crossing.lat <- 88
lat.radians <- star.crossing.lat * ((2 * pi) / 360)
star.crossing.y <- star.radius * sin(lat.radians)
step.length <- .001
planet.x.pos <- seq((-2 * star.radius), (2 * star.radius), step.length) # Position of center of planet position, in terms of stellar radius.
steps <- as.integer(length(planet.x.pos))
point.count <- 25000
box.point.count <- ceiling(point.count * (4 / pi) * 1.1) # Can be wasteful
light.curve <- vector(length = steps, mode = "numeric")
counter <- 1
for(i in planet.x.pos){
#Generate poins on planet
box.points_df <- data.frame(runif(box.point.count, min = -planet.radius, max = planet.radius), runif(box.point.count, min = -planet.radius, max = planet.radius))
colnames(box.points_df) <- c('x', 'y')
box.points_df$r <- sqrt(box.points_df$x ^ 2 + box.points_df$y ^ 2)
box.points_df$in.circle <- box.points_df$r <= planet.radius
circle.points_df <- box.points_df[box.points_df$in.circle == TRUE, ]
circle.points_df <- circle.points_df[1:point.count, ]
#Move planet points into position
circle.points_df$x <- circle.points_df$x + i
circle.points_df$y <- circle.points_df$y + star.crossing.y
#Check if points intersect with star
planet.point.check <- sqrt((circle.points_df$x) ^ 2 + (circle.points_df$y) ^ 2)
intersect.points <- sum(planet.point.check < star.radius)
#Calculate relative flux
pct.planet.intersect <- intersect.points / nrow(circle.points_df)
area.intersect <- pct.planet.intersect * pi * planet.radius ^ 2
relative.flux <- ((pi * star.radius ^ 2) - area.intersect) / (pi * star.radius ^ 2)
light.curve[counter] <- relative.flux
counter <- counter + 1
}
And here is a sample of the output for a transit for the parameters in the code above.
